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JNTRODUCTION 


Let me ask the question most frequentlv asked of owners of 
home computers: “What on earth do you do with it?" Recently, mv 
answer has become: "I use it as a learning toy." This book is 
based on the premise that the home computer is best used as a 
learning device; and 1Е уоц'ге пот careful, vou'll find yourself 
spending hours with it, learnirg again all of those things that 
you had forgotten from school years past. The reader and user of 
this book will experiment with programs and subroutines that inte- 
grate functions, invert matrices, find roots of functions, plot 
curves and contours of surfaces, and solve differential equations. 
He mav find renewed interest and greater understanding of the 
mathematics behind some of the widely-used computer algorithms. 

He тау learn some new programming techniques, since many demon- 
stration programs are followed Бу suggestions for altering the 
main program and/or subroutines to produce customized results. 
And he mav even find some genuine uses for some of these programs 
and subroutines, (as have І). 

The book is based upon 25 programs and subroutines, written 
іп BASIC for the Sinclair ZX81, апа TS2068 computers. The programs 
and subroutines are presented in 59 demonstration problems. No matn- 
ematical fluency or ability is required on the part of the reader 
to try out any of the demonstration programs: each is presented in 
a step-by-step table of operations. 

Many of the central programs and subroutines call upon other 
subroutines within the central set of 25 programs and subroutines. 
In addition, the main programs of each demonstration problem consist 
mainly of subroutine calls to central subroutines. Іп this manner, 
much repetition ከ85 been spared. For example, ROOTN, the root- 
finding program found in chapter 6 calls MXINV, the matrix-inversion 
subroutine found in chapter 5. 

The central programs and subroutines have compatible line- 
numbers, so that they all may be loaded into memory at the same 
time. The first subroutine, VINAX, a subroutine to find the min- 
imum and maximum of a data-vector, starts at line-number 1000. All 
of the demonstration main programs are found between line-numbers 
1 and 999. Тһе central programs and subroutines are numbered 2 
apart, so that modifications to programs may be easily made. 

I have attempted to make the programs and subroutines as gen- 
eral as possible. For example, the linear least-squares regression 
subroutine, REGRP, found in chapter 6, is designed for any general 
sum of fitting functions, (powers of x, exponentials, etc.). Lin- 
ear regression is a sub-category of these procedures and it is pre- 
sented as a demonstration program only. 


The book is organized as follows: Each chapter starts with а 
short and corny introduction just to start things off. Next, the 
program/subroutine listings are given, followed by a "spec-sheet" 
giving the necessary inputs to the subroutine, an example call, and 


a description of the algorithm, ог а line-by-line description of 
the subroutine or program. Several problems are then presented in 
this form: A statement of the problem, a listing of a main program 
to enter, a table of subroutines and programs to enter, а step-by- 
step set of instructions to use the demonstration program, a dis- 
cussion section to give the user some feedback as to whether his 
main program is working or not, and finally a suggestion section to 
give the user some motivation to modify the programs or subroutines, 
or invent new algorithms. 

The last step in each of the demonstration programs is the 
theme of the book: to get the reader/user involved in programming 
and problem-solving by modifying the programs to his/her need or 
desire. There are a large number of infinitely more sophisticated 
programs which run faster or more accurately for specific problems 
such as those in the demonstration programs. However, the programs 
and subroutines in this book are a good collection for general pur- 
poses. The reader who needs a faster or a more accurate approach 
should turn, (quickly), to the reference section. In addition, I 
believe there is enough documentation with each program and sub- 
routine to enable the serious user to modify the routines to more 
specific purposes. For example, the general Simpson's Rule Inte- 
gration subroutine, INTEG.1, found in chapter 4, might be modified 
so that the subroutine returns with an error if the computed 
estimated error gets larger than a given amount. 

Several appendices are included at the end of the book: One- 
line functions, Physical and Mathematical constants, Units con- 
versions, and a table of central programs and subroutines. 


I own a Sinclair ZX80 computer, upgraded with an 8K ROM, and 16K RAM. 
I also own a TS2068. I acquired both of these before TIMEX dropped 
out of the home-computer market, but regret neither purchase. The 
ZX80 doesn't have а "SLOW" mode for flicker-free graphics. Conse- 
quently, none of the programs in this book make use of the "SLOW" 
mode found on the 2Х81. In addition, the changes needed to get the 
programs to run on the TS2068 are minor ones. The color-graphics 
and other enhancements available on the TS2068 are not used here, 
in an effort to keep the programs and subroutines as compatible as 
possible between the ZX81 апа TS2068. 

It is possible to load into 16K of memory all of the programs 
and subroutines in this book, with enough memory left over to write 
a main program. This is the most time saving method of using the 
book: each demonstration problem requires the user to enter at 
least one of the programs or subroutines, and if they are already 
in memory, that much time has been spared. In addition, any of the 
other routines may be used in modified main programs. 

The Sinclair computers offer a'superior BASIC and monitor, in 
my opinion; if not in speed, then in style. It is possible to 
dimension an array by a variable: DIM X(N) is acceptible. This 
book takes advantage of the VAL function, which evaluates string 
expressions to a floating-point number. This is used here for par- 
ameter passing, or in place of defined functions, (which the TS2068 
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has, FN, but the ZX81 does not). Бог example, the subroutine 
MXINV inverts the matrix named in the string Xs. If X$="ACI,J)", 
then MXINV will invert the matrix АС mE 

The monitor program provides Single keyword entry. I have 
added a full-size keyboard to my ZX80, and have found that, with 
practice, it is possible to touch-type these Keywords in, saving 
a lot of time in program entry. Two more very attractive features 
of Sinclair BASIC are its syntax-checking and the ability to 
CONTINUE а program even after editing program lines. 

The central programs and subroutines can be easily translated 
to any other BASIC environment, with minor changes. I have trans- 
lated the matrix-inversion subroutine МХІМУ to a Commodore VIC20 
with no difficulty. However, some things to watch out for when 
translating: Оп the 2Х81 and TS2068, DIMensioning an array sets 
all values to 0; VAL on other machines tend to accept only string 
expressions involving one floating point number only. | 


Beginning programmers might find the chapter on plotting the 
most interesting one to start out with. Follow the instructions in 
the demonstrations, and then experiment with уоцг own variations on 
theme. Move оп to the last chapter апа try out some multiple- 
precision additions, subtractions, and multiplications. (Try your 
hand at writing a routine for multiple-precision division). 

High-school/college students might use the book as a Workbook: 
to get a better picture of what.he is trying to learn in Algebra, 
Vector-analysis, and Calculus. Use INTEG.2 to numerically inte- 


tne 


grate those impossible homework Problems, (only as a check, of course.) 


The person who really needs to find roots to equations, find 
best-fit curves, or solve differential-equations may find the 
programs in chapters 6 (Roots and Regression), and 7 (Differential- 
equations? adequate, or he/she may want to translate them to 
higher-powered machines. 


But to all users of this book: Have a good time experimenting 
with the programs and subroutines, for it was in that Spirit that 
this book was written. I would like to hear from anyone who has 
improvements/additions/suggestions/additions/translations/(and any- 
thing else.) 
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USING THE TAPE 


The tape supplied with this book holds two copies of the file 
named "DEMOS", which contains all of the demonstration programs апа 
central subroutines and programs for the TS2068. There is enough 
space left on the tape to record your own programs. 

To use the tape, type 


LOAD "DEMOS" or LOAD ”" 


and start your recorder. The program takes about 24 minutes to load. 
Once the program has loaded, it will start automatically. A Title 
page is first displayed for a few seconds, followed by a menu of 

all of the demonstrations, (answer "y" to the "scroll?" Prompt to 

get the entire menu listing). Enter your choice by typing the number 
of the demonstration you want to run. Find the demonstration in the 
book, and proceed as shown in the "Run program" step. 


USING THE EOOK 


Each chapter in the book is presented in this format: 
1. Introduction 
2. Central subroutine(s) or program(s): 

a) Listing(s) 

Ы) "Spec-sheet(s)" 
3. Demonstrations 


The "spec-sheet" lists the properties of a particular subroutine or 
program. It has this format: 
1. Short description of what the subroutine or program does. 
2. A flow-chart of subroutine or program use: 
a? For a program: 
RUN: 
(program name): List of prompts, and variables changed. 
STOP: Expected results. 
b? For a subroutine: 
CALL: List of variables which must be initialized or 
dimensioned before entering subroutine. 
(subroutine пате): List of prompts, and variables changed. 
RETURN: Expected results. 
3. Description of algorithm or line-by-line description of program 
or subroutine. 


The demonstrations illustrate the properties of the central 
subroutines and programs. To use a demonstration, follow the step- 
by-step instructions to enter the main program, central subroutines, 
and run the program. The step for running the program has two columns. 
The idea is to enter the USER ENTERS item followed by an {enter}, and 
wait for the COMPUTER RESPONDS item. Use the "discussion" and 
"suggestions" steps for your own programmingexperiments. 


VECTOR OPERATIONS 


Chapter 1 


Just аз some of us do and some of us don't have a direction 
in life, vectors do and scalars don't have an associated direction. 

For the time being, let's assume we are talking about 2-dim- 
ensional vectors. Any general (2-d) vector may be described by 
two "coordinates" x, апа x,, the x- and y-components of the vector 
with respect to somé referénce system, as shown in figure 1-1. 
In figure 1-1, the vector is "associated" with the point (a,b). 
The vector might represent the instantaneous velocity (speed and 
direction) of a neutron, at a particular time when the neutron 
is passing through the point (a,b). 

The direction of Х in the figure seems to be about NE, (if N 
is top). X may also be described by two other numbers: its length 
and its angle with respect to the +X-axis, where 


7—5 = 


LENGTH(X) = хі “*> 
ANGLE(R) = ATN(X,/X,) 


Moving up to higher dimensions, it is clear that a set of three 
coordinates, x „Хо, Ха» describes a certain vector іп three dim- 
ensions. Its tength and angles with respect to a given coordinate 
system may also be defined similarly to the 2-d case. But there 
are sets of 4,5,...(even о) numbers which may also be regarded аз 
vectors in dimensions which can only be imagined. 

This chapter presents two subroutines which manipulate the 
coordinates of vectors. As in the case of the 2-dimensional vector, 
lengths and directions are determined by use of the subroutines. 

The subroutine VINAx finds the minimum and the maximum of a 
vector's coordinates. This is used in the next chapter by scaling 
subroutines, which are, in turn, used by plotting subroutines. 

The subroutine VDOTP can be used to perform various summation 
operations. The dot-product between two vectors is one such opera- 
tion. 

Four demonstration problems are given here to demonstrate the 
vector operations subroutines. In larger programs, it may be 
advantageous to consolidate the subroutines into the main program. 
For example, the polynomial regression subroutine, REGRP, found in 
chapter 6, could have used many calls to the vector operations sub- 


routines; these have been consolidated into one main summation loop. 


LENGTH(X) 


ж 
А, 


ANGLER) 
І 


ae figure 1-1 Vector 


јада БЕН VINA ПРЕСТОВ WIN, MAXI 
1002 LET Юзі1Е37 

1004 LET ኒ|=>፦ህ 

1982 FOR 151 TO hi 

1088 LET Hz'RL = 

1210 IF НА) THEM LET ЕМ 

1012 IF Н:Ч THEN ШЕТ ЏЕМ 


2814 па Т 
1188 BEN UDOTP (VECTOR DOT-PRODU 


1198 NEXT. 


Figure 1-2 Vector Operations Subroutines 


УІМАХ 


This subroutine returns the minimum апа maximum element of а 
vector named by X$. 


CALL: 


VINAX: 


RETURN: 


X$ = name of vector, indexed by I 
N = Dimension of named vector 


Changes 1 М 0 М 


U = Maximum element of named vector 
V Minimum element of named vector 


Example call: 


LET VINAX-1000 
LET N-10 
DIM УСМ) 


LET X$="Y(1)" 
GOSUB VINAX 


VDOTP 


This subroutine returns the dot-product of two named vectors. 
It may also be used to compute other types of summations, such as 
the sum of the elements of a vector. 


CALL: 


VDOTP: 


RETURN: 


X$ = string naming summation operation with vectors 
indexed by I 

N = Dimension of named vectors 

Changes I V 


У = Dot-product (or summation) 


Example call: 


LET VDOTP=1100 
LET N=10 
DIM X(N) 
DIM Y(N) 


LET X$s"X(CIOxY(IO" 
GOSUB VDOTP 


є 


«ሀ fpe P e p e p p E E ра ЛЕ» 


Demonstration 1.1 VECTOR NORMS - 


Step 1 (PROBLEM) : Print the values of 


these vector norms of a 


vector X=(X(1), .,XCN)): 

id 
Euclidean norm: Xile = X5۶ 
Maximum norm: |[|ቿ|| = тах(|Х(1›|) 
Summation norm: | ||ቿ||ፊ = Eos 

~ ај 


Standard viation: б- 


for the vector X = (-1,2,5,3,-1) 


Step 2: Enter main program: 


180 >БЕМ DEMOl1:1 
101 LET VINAx=laae 
162 LET ህይርፐጅ =3386 
185 LET ከ=፳ 


z I 
UDOTE ы 
PRINT "EUCLIDEAN NORM = "IR 


hire rà 
ка 

лс 
ሂስ 
С? 
га 


© г; 


ሁሁ вена 


| Thum NORM = "U 
PRINT ра = “;ህ 


жи 
S d 
са 
ея 
и пи 
- _ ዘ 
25x 
КЛ —. 
ኀ 
a — 


KES (X TIME ee CD) -ቨ1” 
SUB ЧПОТЕ 

Т "MEAN = 7:84 

T "STANDARD CEVIATION = 


"Mi 


“оо 
гей F C9 Р Fam eo ^J Cn en p ra 


пи руга 


[a] 
iino на Г в በር) 


пу- 


tep 3: Enter subroutines:  VINAX VDOTP 


Step 4: Run program: 


USER ENTERS: 


R 


N 


СХ СІ): -M)* ; where mean: M 


COMPUTER RESPONDS: 
(results printed) 


Step 5: Discussion апа suggestions: 


The norm of a vector can be thought of as its "length", or at 
least a characteristic scalar value associated with the vector, by 
which it may be compared with other vectors. 

Any norm must obey several criteria: 

СИ | is the norm of the vector Ж) 


1. ||ቿ|| >= о; 1111 = O if апа only if X = Û 
2. |I*-Y11 <= [IRI] + LIVI 
з. llaX|| = 1|8|።!|!|ጆ!|| 


The four norms іп this demonstration each assigned different 
numbers to the vector Я-(-1,2,5,3,-2), but each still obey the above 
criteria. 

Be sure to try some other vectors besides the particular one 
given in the demonstration. Try vectors of different dimension 
also, by changing lines 103-104 to the appropriate dimension). 

In particular, try the zero vector and thereby verify one of the 
criteria for vector norms. In fact, try to verify all of the other 
criteria, by examples. 

One problem with this demonstration is that it is not general: 
the vector dimension and vector coordinates are not changeable 
without modifying program lines. In the next chapter, a subroutine 
is presented for entering data points. Using this idea, write a 
more general demonstration program. 

Find some other vector norms, or invent some, and verify the 
criteria for vector norms by examples. 


Demonstration 1.2 SCALAR OPERATIONS BETWEEN TWO VECTORS 


Step 1 (PROBLEM): Print the values of these operations between 


vectors X-(X(12,...,XXOD) and ў= (У(1),...,УСМ)): 


ч 
Euclidean dot-product: (፪,ኛ) = T ad 


Covariance: соу(3,%) “Т Cac -My ye СУ СІ)- zM )) 
ң T54 


where: M 1/8) ፳(፤) 
x ri 


M = 1/N YOD 
y 2-4 
Distance between Х ара У: N 
баіз%(5,%) = ) aon УСІ)» 
Апа1 etween Х апа Ү (іп degrees): 
angle(X,9) = ACS £X.) х180/9 


КА, ፍ(ነ. 


for the vectors X з (-1,2,5,3,-1) and Y = (2,5,7.-1,40) 


Step 2: Enter main program: 


138 REM РЕНОЈ.2 
121 LET UDüOTPziicg 
132 LET HzE 

133 DIM RIS! 

124 DIM У (5) 

135 LET жідіст-і 
138 LET XxíZ!sZ 

137 LET KESI -5 

i38 LET хійзш3 

138 LET Х#[51 2-1 
142 L “дъвка 

141 LET Wiis 

142 LET Vidi ат 

143 LET YI] 2-1 
144 LET WiSisda 
145 LET Хфа"х (Liev (І) 
146 GO SUB ПОСТЕ 


Iso 52 sib vboTe 
151 LET HX=UzN _ 


153 


ÉET ዘሃ=ህ።ክ 
LET X&z" ХІІІ -МХІ ЗСУ ITI HY 


к 
ле 


GO SUE uDOTFE 
PRINT "СОМІХ,Ч) = 
LET Хай туст =ፐ el) 


GO SUE ህርዉፐጾሯ 

PRINT "РІЗФТІХ, У)! = "58 У 
LET Ж ж CT жету 

GO SUE ОТР 

LET Xxx =u 

LET Хаш“ҮШІржҮ(СІ!" 

ca SUE чосте 
LET ኙሃ =ህ 
PRINT "ANGLE ( 
ወጾ хаттау 


> STOP 


+, Са Сид Cri Ch Ti Chi С СП mid eai 


нашр чи uia 


“PEEP “ሀ ። s HF: HE 


їж 
p 


- 
т 
th 


Step 3: Enter subroutines: VDOTP 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN (results printed) 


Step 5: Discussion and suggestions: 


The four operations in this demonstration each assign a 
different number бо the two vectors X and Y. There are two 
"geometrical" operations:  dist(X,Y32, and angle(X,Y5. For two- 
dimensional vectors, the geometrical picture is this: 


Try different vectors and dimensions. Be sure to try some 2- 
апа 3-dimensional vectors to look at the "geometrical" operations 
for "real-world" vectors. 

What do the operations become for 1-dimensional vectors? 
For example, the Euclidean dot-product degenerates to multiplication 
for vectors with only one coordinate. | 

Find some other operations between vectors, ог invent some: 
there are an infinite number of possibilities). All of these oper- 
ations give a scalar result (a real number, as opposed to an 
ordered set of real numbers, that is, a vector). There are also 
operations, such as vector addition, which give a vector result. 
Write your own subroutines to take care of any operations between 
two vectors which aren't covered by VINAX or VDOTP. 

Finally, the demonstration's main program might be made more 
general by the addition of a portion for entering dimensions and 
coordinates of vectors. A routine similar to the subroutine 
POINT.1, presented in chapter 2, could be added in the place of 
lines 132 to 144 of the main program. 


11 
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Demonstration 1.3 DISCRETE CONVOLUTION 


Step 1 (PROBLEM): Plot the discrete convolution of Х and У: 


Ү(І). = (о for І<-0 
1 for 0<1<=60 


(Y is a unit step function.) 
KCI) = 0 for І<-0 
20/1 for 0<1<=60 


Step 2: Enter main program: 


-4-4-4-4-4-4 


каб у ПИ FE MME 


= 
id 
=: 
г 
= 
E 
= 
Ё 


Step 3: Enter subroutines:  VDOTP 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN (wait about 55s) (results plotted) 


Step 5: Discussion: 


The discrete convolution given in this demonstration is the 
discrete version of analog convolution. Analog convolution can 
be described as follows, (the example corresponding to the 
demonstration). 

Suppose some system gives an output signal when it is presented 
with an input signal. The "system" might be an electronic circuit 
and the "signals" voltage waveforms. If the system is presented 
with a unit impulse (a voltage spike in time, with 


со 
fat Yay =4 
that is, unit area), it responds with the output function 


X(t) = 20/(t+1) 


20 
* X 
Y 
system 
Е 0 t 0 


The function X(t) is called tne system's impulse response. ТЕ 
the system із presented with another sort of input: 


Y(t) = unit step function 


the system's response will be the step-response of the system. In 
general, the output of the system will be the result of convolving 
the input signal with the system's impulse response: 


з Op 
Convolution(X,Y) = [xcsrvce-sras 


«є 
So the step-response is given Бу: 
t 


ЕЕЕ = 20*LN(t+1) 
0 
У - - X 
system 
t О t 


Step 6: Suggestions 


It may be possible to generalize this program to accept 
functions in the form of strings. See chapter 6 for some hints. 

The function X(I) may be varied by changing the constant 20 
in line 176 of the main program. See what effect this has on the 
results. Try other impulse response functions, Х(1), and system 
input functions, УСІ). 


18 


‘Demonstration 1.4: GRAM-SCHMIDT ORTHONORMAL IZATION 


Step 1 (PROBLEM): Е1па an orthonormal basis for the vectors 


X,, Ху, and X, ек 
x, вс, Е POND) 
x, = t0 TLD 
X, = (0,0,1,1,1) 


Step 2: Enter main program: 


3 
z 
E 
z 
5 
5 
5 
а 
8 


fan eu ЕЛ РУКА EO IP ба D Cü 4 On (በቶ ra D P O 


I2 ES EIE ES ВО E Rs ра Ea FI јајета 


ue» 
I p SIG O СТС E La fO I 


1 
1 
1 
1 
፲ 
1 
1 
1 
1 
1 
1 
5 
= 
= 
= 
> 
Е 
2 
= 
= 
= 
= 
= 
= 
2 
Е 
= 
= 
= 
= 
= 
= 
= 
= 
E 
z 
= 
Е 
2 
Е 
2 
= 
= 
1 
= 
= 
= 
= 
= 
= 
5 
= 


Step 3: Enter subroutines: VDOTP 


Step 4: Run program: 


USER ENTER: COMPUTER RESPONDS: 
RUN (wait about 10s) (results printed) 


Step 5: Discussion: 


This demonstration finds а set of basis vectors (Е), 5- 
dimensional vectors, which have the following properties: 


i. Each Ё has a Euclidean norm of 1 (see demonstration 1.1). 
2. The Euclidean dot-product (Ё.,Ё »'= 1 if ፲=ጋ ог О if І‹Ј. 
This means that the basis vectors are orthogonal to one another, 
the angle between any two is 90 dggrees (see demonstration 1.2). 
3. Any vector in the vector-space ІВ", the collection of all 
5-dimensional vectors, which can be written as a linear sum of 
X- d Хо. апа X, (given in the demonstration), may also Бе 
written as a linear sum of the basis vectors {E,). That is, if 
ኛ can be written аз the sum: 
=. =~ - -ь 
= a,X, + وگو‎ + Яка 


Then there also exist b's such that Y can also Бе written: 


Ү = Е: + Е. + bE, 


Note that, in the print-out section of the main program, lines 
215 to 221, the results are truncated (line 218). There still 
may be trouble for some other vectors if the truncation error is 
less than the actual coordinate value. 

In addition, the Gram-Schmidt method requires that the vectors 
(Х,) be linearly-independent. That is, no one of the X's пау be 


written as the linear sum of the rest of the R's. 


Step 6: Suggestions: 


Lines 192 to 199 might be replaced with a routine for entering 
any set of vectors from the keyboard. 

Analyze the workings of the program with a flowchart, and note 
any difficulties which might occur with the algorithm for certain 
cases. For example, what if the number of vectors is greater 
than the dimension of the vector-space? This set, of course, would 
be a linearly-dependent one. If possible, correct some of tne 
deficiencies in the program. 

as always, try other vectors than those given in the demon- 
stration. In particular, let one of the X's be'the zero vector. 
(The zero vector сап always be written as the linear sum of any 
other set of vectors). 
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- Chapter 2 PLOTTING 


There comes a time in everyone's life when he or she must 
deal with data. The person who wishes to have a picture of what 
his data is doing plots it. It is to this person that this chapter 
is addressed. Іп addition, any person who finds watching patterns 
a worthwhile activity should listen in. 

The subroutine PLOTD plots data points, which are input via the 
subroutines POINT.1 or POINT.2.  PLOTD also makes use of SCALE.1 or 
SCALE.2 to scale the data vectors so that all of the data points 
fit within the plot.  PLOTD assumes that the data it is given is 
cartesian in nature, (х- and y-coordinate). Once it is called, 
it finds an appropriate plotting window via SCALE.1, or SCALE.2, 
and then plots the data found in vectors XO and YC) in the form 
(ХСТ),УСІ)). The difference between SCALE.1 апа SCALE.2 is that 
SCALE.1 finds endpoints which may be expressible as an integer times 
a power of ten, where SCALE.2 finds the minimum and maximum of the 
data vector, (ХО, or ҮС»). 

Data is supplied to PLOTD via one of two subroutines, POINT.1 
and POINT.2. POINT.1 prompts the user (you) for the number of 
points and each data point, (in which case you may either respond 
or go back to sleep.)  POINT.2 prompts the user for an expression 
of X for sampling at equal intervals. POINT.2 also asks for how 
many points to be gathered, and over what range of X values to 


sample. 
Demonstrations 2.1, 2.2, and 2.3 illustrate the uses of these 
subroutines. In addition, demonstrations in other chapters make 


use of some of these subroutines for their particular purposes. 
The demonstrations in this chapter plot cartesian data, polar data 
and a sampled function. 

The program PLOT2 is a program that plots plane-curves which 
are described by any of 4 different types of curve descriptions: 
cartesian, cartesian parametric, polar, and polar parametric. 

The user supplies PLOT2 with an expression or expressions to be used 
to generate a plot, the desired plotting window, and for parametric 
and polar plots, a range for the parameter to be swept over. 

The program PLOT3 generates a picture of tbe contours of a 
three-dimensional surface. Тһе user supplies en expression іп 
X and Y, a plotting window, and a parameter called the contour- 
spacing-parameter, (C-S-P). The C-S-P determines how far apart 
to space the contours. 

Example curves and contours are given, for comparison with 
results from the demonstrations. There is a lot of room for experi- 
mentation in this chapter. Changing the value of a constant or two 
in the expressions of the demonstrations can give very different 
results. 

Plotting usually takes considerable amounts of time, and you 
may want to do it in conjunction with another activity, such as 
writing a letter to your congressman. For example, it took about 


13 minutes to plot the contours of the magnitude of the complex 
sine function, (on the 2X81). Not all of the plots in the demon- 
strations in this chapter take quite that long. However, be patient! 
Using a computer to do your plotting is still much faster and easier 


than doing it by hand. 
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Figure 2-1 Data point input, bounding, and plotting subroutines 


POINT. 1 


This subroutine prompts the user for data points. The number 
of data points is first entered, then each point is entered into 
vectors ХО) and YO). 


CALL: 


POINT.1 Changes І М XD ҮС) 
Prompts, in order of appearance: 
1. "INPUT NO. PTS." 
2 "СТО SC ys. № 


ХО - x-coordinates of each data point 
YO) + y-coordinates of each data point 
М = dimension of X() and У() = number of data points 


RETURN: 


Example call: 


LET POINT-2000 
GOSUB POINT 


Summary of program operation: 


Lines 2002 to 2008: prompt user for number of points, and dimension 
arrays X(N) and У(М). 

Lines 2010 to 2022: prompt user for each point and enter values 
into Х(1› апа уст», 
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POINT.2 


This subroutine prompts the. user for an expression to ђе sampled 
at equally-spaced intervals. The expression is entered in the form 
of a String variable, as a function of X; for example Y(X)-"Xx*x2-«3/X". 
The number of data points is first entered, then the expression is 
entered, and then the range of X-values for sampling is entered. 
The subroutine then samples the function and enters the x- and y-coord 
inates into the vectors ХО) and YO). 


CALL: 


POINT.2 Changes A B М XO YO Y$ 
Prompts, in order of appearance: 
1. "INPUT NO. PTS." 
2. "INPUT Y EXPRESSION-Y(X)O" 
3. "INPUT XMIN,XMAX" 


RETURN: ХО = x-coordinates of each data point 
Ү() = y-coordinates of each data point 
dimension of ХО and У() = number of data points 
- Xmin 
B Xmax 
УФ = Y expression, Y(X) 


> 2 
ዘ 


Example call: 


LET POINT-2100 
GOSUB POINT 


Summary of operation: 


Lines 2102 to 2108: prompt user for number of data points, and 
dimension X(N) and Y(ND. 

Lines 2110 to 2118: prompt user for a function of X, and Xmin,Xmax. 
Lines 2120 to 2128: sample the function and enter the samples into 
ХСІ) апа УСІ»). 


SCALE.1 


This subroutine finds lower and upper bounds for a vector 
named іп X$. The lower and upper bounds are each an integer times 


a power of ten. 


CALL: VINAX = 1000 
X$ = Name of vector to be bounded, indexed by I 
М = Dimension of vector 


SCALE.1: Changes I M PU V 
Calls VINAX 


Lower bound 
Upper bound 


RETURN: v 
U 


Example call: 


LET VINAX=1000 
LET SCALE=2200 
LET М-10 
DIM X(N? 


LET X$="X(I)" 
GOSUB SCALE 


Summary of operation: 


Line 2202: determine minimum and maximum values of data vector 
named by string X$. For example, if X$="K(I)" then the minimum 

and maximum of X() are found. 

Lines 2204 to 2212 provide for the case that the minimum and maximum 
of the data vector are the same. In this case, the returned 

lower and upper bounds are separated by 1. 

Lines 2214 to 2228 determine the lower endpoint. 

Lines 2230 to 2238 determine the upper endpoint. 
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SCALE. 2 


This subroutine finds lower and upper bounds for a vector 
named in X$. The lower and upper bounds are the minimum and 
maximum values of the data vector. 


CALL: VINAX = 1000 
Х$ = Name of vector to be bounded, indexed by I. 
М - Dimension of vector 


SCALE.2: Changes I MU V 
Calls VINAX 


Lower bound 
Upper bound 


RETURN: V 
U 


И 


Example call: 


LET VINAX=1000 
LET SCALE=2300 
LET N=10 
DIM X(N) 


LET X$-"X(I)" 
GOSUB SCALE 


Summary of operation: 


Line 2302: Determine minimum and maximum values of the data vector 
named by string X$. 

Lines 2304 to 2312 provide for the case that all values of the 
data vector are the same. The returned lower and upper bounds, in 
this case, are separated by 1. 


PLOTD 


This subroutine plots data in the vectors XO and YO. ХО 
contains the x-coordinates апа У() contains the y-coordinates of 
each data point. For example, data point number 5 is 


X = Х(5) ; ሃ Y(55 


2200 for SCALE.1 or 2300 for SCALE.2 
1000 

ХО - x-coordinates of each data point 
y-coordinates of each data point 

М = Number of data points 


CALL: SCALE 


« 

= 

= 

> 
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PLOTD: Changes А В C ር І M N P U V XO X$ YO 
Calls SCALE.1 ог SCALE.2 (which call VINAX) 


RETURN: Data plotted 
Lower left and upper right endpoints of plotting-window 
printed at top of screen 


Example call: 


LET PLOTD=2400 
LET VINAX=1000 
LET POINT=2000 
LET SCALE=2200 
GOSUB POINT 
GOSUB PLOTD 


Summary of operation: 


Lines 2402 to 2416: Determine plotting window. 

Lines 2418 to 2424: Plot each data point. 

Line 2426: print lower left and upper right end-points of the 
plotting window. 


23 


Demonstration 2.1 PLOT CARTESIAN DATA 


Step 1: Enter main program: 


Jm p ram у 


Step 2: Enter subroutines:  VINAX POINT.1 SCALE.1 PLOTD 


Step 3 (PROBLEM): Plot the cartesian data in the table below: 


1: кс): мета: 
1 D о 

2 -1 -1.6 
3 -.36 -.4 
А „2 „14 
5 Зи 25 
6 1.7 38 
7 2.1 52 
а 3.5 O 


Step 4: Run program: 


USER ЕМТЕВЗ: COMPUTER RESPONDS: 


RUN INPUT NO. PTS. 


8 X(1),Y(1)=? 
S. u uuu. ی ی ت چت ت‎ uU ማመ 


0 Х(2),У(2) =? 


enter дата from table in step 3 


2 K(8),Y(8) =? 


(results plotted; 


24 


Step 5: Discussion of results 


The data in this example was fabricated. Any resemblance to 
actual data is purely coincidental. 

Line 231 of the main program selected POINT.1 to fill the 
data vectors X() and Y(), so that vou could hand enter the data. 
SCALE.1 found two endpoints for each data vector so that all of 
the data would appear on the plot. Since SCALE.1 was used, the 
two endpoints for each axis should be integers times а power of 10. 
In this case: 


Xmin = -1 ; ( = -1ж100) 
Хтах з А 

Ymin = -2 

Ymax = 1 


The lower left and upper right endpoints of the plot are 
printed above the plot. In this case: 


Lower left: (Xmin,Ymin?) = (-1,-2) 
Upper right: (Xmax, Ymax) = (4,1) 


PLOTTING WINDOW: 


x} (Xmax, Ymax) 


(Xmin,Ymin) |х 
Step 6: Suggestions 


Since you have spent some time entering the subroutines and 
main program, SAVE them onto a file for future use. You may want 
to keep the subroutines in memory for the time being and try demon- 
strations 2.2 and 2.3. 

This demonstration can be used to plot any other cartesian 
data you may have lying around. Hint: look in a magazine for an 
article on the economy. Chapter 6 also contains some tables 
of data which are used in the demonstrations there. If you don't 
have any data on hand, make some up and plot it. Generate some 
figures in this way. 

For a variation on this demonstration, try to write а main 
program which makes use of PLOTD, or а similar subroutine, and which 
plots one data point at a time. The plotting window would have to 
be pre-determined, and SCALE.1 and SCALE.2 subroutines wouldn't be 
needed. Since PLOTD calls SCALE automatically, you can get around. 
this by writing a subroutine SCALE.3: 
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2700 
2702 
2704 
2706 
2708 
2710 
2712 


REM SCALE.3 

LET U={pre-determined value} 
LET V={pre-determined value} 
ЈЕ Хф-"Х(Т)" THEN RETURN 
LET U={pre-determined value) 
LET V-(pre-determined value} 
RETURN 


Or more simply, edit out some of the lines in PLOTD. 
For another variation, try to combine the main program and 
subroutines all into one big program. 


Demonstration 2.2 PLOT POLAR DATA 


Step i: 


s 


448. 44443 


arm РРА rm t 


"rraworrr rm 


RETI Fr CH 


Step 2: 


г 
IU 


орні я 
M 


Enter main program: 


Enter subroutines:  VINAX POINT.1 SCALE.2 PLOTD 


Step 3 (PROBLEM): Plot the polar data in the table below: 


155 хе! RCI 
1 O° о 
2 15° 52 
3 45° 26 
4 50° 25 
5 70? а 
6 75° 21 
7 80° ዕ 
Step 4: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT NO. PTS. 
7 X(1),Y(1)=? 
0 
0 Х(27,У(2)-? 
enter data from table іп step 3 
£i - Х(7),Ү(7)=? 
80 
ዐ undi (results plotted) 
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Step 5: Discussion: 


You should have obtained a plot similar to figure 2-2, below. 
This is a plot of the points in the table on polar graph paper. 

The demonstration is an example where, after data vectors X() 
and Y() have been filled via POINT.1, the data is manipulated 
(lines 246 to 250) before it is submitted to PLOTD for plotting. 
Thus, in this example, the polar data entered via POINT.1 is 
converted to cartesian data, which is the appropriate form for 
PLOTD. Note that in line 247, a degree to radian conversion 
is made, (see Appendix А). 

Line 243 selected SCALE.2 for finding appropriate endpoints 
for the X and Y axes. The endpoints should be the minima and maxima 
of the data vectors ХО) and 70. 


Figure 2-2 Polar plot 


Step 6: Suggestions: 


If you have need of а polar data plotting program from time 
to time, SAVE the main program and subroutines in this demonstration. 
Before you throw away what vou have done with a NEW, try plotting 
some different polar data. Create some data of your own or use а 
function from table 2-1 to generate a few points. You may want to 
keep the subroutines in memory for the next demonstration. 

For a variation of this example, alter the main program so 
that it plots polar data in the form of radians and radii. 

You can insert another type of data manipulation routine 
before plotting the results. For example, a derivative-taking sub- 
routine, DERIV.1 is found in chapter 4. The demonstration of DERIV.1, 
demonstration 4.2 differentiates the data before plotting. 

For another variation, combine the program and some or all of 
the subroutines into one main program. 


Demonstration 2.3 PLOT SAMPLED FUNCTION AS POLAR DATA 


Step 1: Enter main program: 


те 


LOTT ከበ9)][”፦ታ 


й ура 


ርበ 
+ 


+ ГА ГА ГА EF] C3 D ጣጠጠጠጠበ 


Әсер 2: Enter subroutines: VINAX POINT.2 SCALE.2 PLOTD 
Step 3 (PROBLEM): Sample the function below and plot as polar data: 
Есе» = СО5(8) 3 Ос <а< + 41 


Step 4: вип ргодгат: 


USER ENTERS: COMPUTER RESPONDS: 
ج‎ PRO 
200 INPUT Y EXPRESSION=Y(X) 

COS X я INPUT ХМІМ,ХМАХ 
О 
РІ (wait about 1m) (results plotted) 


Step 5: Discussion and suggestions: 


You should have obtained a plot of a unit circle with center at 
С сб ДЫЗ m 

The program PLOT2, listed іп figure 2-3, сап perform the same 
function as this demonstration program: plotting a polar function. 
However, the use of PLOTD does have the advantage in that after a 
function Пав been sampled, the data may be manipulated before sub- 
mitting it to be plotted. 

In this example, POINT.2 has been used to sample the function at 
equal intervals, and SCALE.2 has been used to give the maximum plotting 
resolution (SCALE.1 would have found a somewhat larger plotting | 
window). 

You will find more examples of polar functions to sample in 
table 2-1. Try plotting one or two of them. 

For a variation on the theme of this demonstration, look at 
demonstration 4.2, where some data has been differentiated before. 
sending it off to be plotted. 


29 


30 


свой REH PROGRAM FLOT2 (PLOT CUR 
JES 


2582 PRINT „СВЕТ, CART. PARA. , POLA 
ROR POLAR PBRR,2' 

2504 PRINT “INPUT 1,2,3,0R à" 
2506 INPUT 2 

2508 PRINT “PLOTTING итнроцт" 
2510 PRINT “INPUT XMIN,XHAX,YF 


мах" 
2512 INPUT ñ 
2514 ТЫРЫТ В 
2515 INPUT C 
2516 INPUT D 
25208 [ЕТ EzH 
25223 LET Е-в 
2524 IF Е=1 THEM Gü то 
=5=5 PRINT "PARAMETER. 
ETT 
2523 PRINT "INPUT THIN ТНЕУ" 
2550 INPUT E 
2532 INPUT Е 
2534 IF Z 
Z535 F 
2538 Е 
UT 
= INPUT x$ 
БЕТМ ` 
FRIR 
ri" 
ІНЕ! 
ELS 
FOR 
LET. 
IF = 
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Figure 2-3 Plane curve plotting program 


PLOT2 


This program plots plane curves described by one of the following 
representations: 


1. Cartesian expression: У = F(X) 
2. Cartesian parametric expressions: X = F, €T) ; Y = Е „(Т) 
3. Polar expression: R = F(T) 
4. Polar parametric expressions: R = Е (T> 3 в = Fa (т? 
RUN: 
PLOT2: Changes A B C D E F M T X X$ Y ሃቁ Z 
Prompts, in order of appearance: 
L "САВ. , CART. PARA. , POLAR, or POLAR PARA.?" 
"INPUT 1,2,3,0R 4” 
2. "PLOTTING WINDOW?" 
"INPUT ХМІМ,ХМАХ,ҮМІМ,ҮМАХ" 
3. "РАБАМЕТЕВ/ТНЕТА RANGE?" 
"INPUT ТМІМ,ТМАХ" 
4. "X/THETA EXPRESSION?" 
"INPUT X(T) OR ТНЕТА(Т)" 
5. "Y/R EXPRESSION?" 
"INPUT Y(X) OR Y(T) OR R(T)" 
STOP: Curve plotted 


Summary of operation: 


Lines 2500 to 2546: Prompt user for type of curve, plotting window, 
independent variable range, and X/THETA and Y/R expressions. 

Line 2550: Compute step-size and limits for sweep of independent 
variable, (X or T). 

Lines 2552 to 2556: Assign values То X and Y for each value of T. 
Lines 2558 to 2564: In the case of polar plots, convert from polar 
coordinates to cartesian coordinates. 

Lines 2566 to 2573: Limit coordinates to plotting window. 

Line 2574: Convert from cartesian coordinates to pixel coordinates, 
and PLOT. 
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Demonstration 2.4 PLOT SEVERAL CURVES FROM TABLE 2-1 


Step 1: Enter program: PLOT2 
Step 2 (PROBLEM 1): Plot the serpentine with the cartesian form 
У = Х/((Х/62%-1/9) 


This is shown in figure 2-4 . 
Use the plotting window: (Xmin,Ymin) = (-10,-10) 


(Xmax,Ymax) = (10,10) 
Step 3: Run program: 
USER ENTERS: 5 COMPUTER RESPONDS: 
GOTO 2500 CART. ‚ САВТ. PARA. , POLAR, ОВ POLAR PARA. ? 


INPUT 1,2,3,OR 4 
1 PLOTTING WINDOW? 
INPUT XMIN,XMAX, YMIN, YMAX 


що» o a вна ው 


10 
-10 
10 Y/R EXPRESSION? 


INPUT Y(X) OR Y(T) OR R(T) 


X/(Xx*xX/36*179) 


(wait about 205) (results plotted) 


Step 4 (PROBLEM 2): Plot the Lissajous pattern with the cartesian 
parametric form 


ЗІМ(СТ) 
SIN(3*T/4) 


X 
Y 


for -20<-Т<-20. 
This is shown in figure 2-5 
Use the plotting window: (Xmin,Ymin) 


እ መቁ.“ 1.2 


(Xmax,Ymax) = (1,10 
Step 5: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
GOTO 2500 CART.,CART.PAR.,POLAR,OR POLAR PARA.? 
INPUT 1,2,3,ОБ 4 
2 PLOTTING WINDOW? 


INPUT ХМІМ,ХМАХ,ҮМІМ,ҮМАХ 
= 
1 
=] 
1 PARAMETER/THETA RANGE? 
INPUT TMIN,TMAX 


-20 


20 Х/ТНЕТА EXPRESSION? 
р Q 

ج ن دت 

SIN T Y/R EXPRESSION? ' 


INPUT У(Х) OR УСТ) OR R(T) 


SIN (.75хТ) 
(wait about 355) ... (results plotted) 


Step 6 (PROBLEM 3): Plot the цітасоп of Pascal with the polar equation 
В = 2xCOS(8) + i 


for -4<-Ө<-4 
Use the plotting window: (Xmin,Ymin) = (-.5,-2) 


(Xmax, Ymax) = (3,2) 
Step 7: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
GOTO 2500 CART. , CART. MARA. , POLAR, OR POLAR PARA. ? 


INPUT 1,2,3,0R 4 
PLOTTING WINDOW? 


3 
INPUT ХМІМ,ХМАХ,ҮМІМ,ҮМАХ 
== 5 . 


-2 
2 PARAMETER/THETA RANGE? 

.INPUT TMIN,TMAX S 
-PI . n 
PI Y/R EXPRESSION? 


INPUT Y(X) OR Y(T) OR R(T) 


2*COS Т+1 
(wait about 45s) (results plotted) 
— р (гезоиіз plo 
Step 8 (PROBLEM 4): Plot the Rhodonea with polar parametric equations 
В = COS(T) 
8 = T/7 
for -%<=T<=4 
Use the plotting window: (Xmin,Ymin) = (-1,-1) 
(Хтах,Үшах) = (1,1) 
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Step 9: вип ргодгат: 


USER ENTERS: COMPUTER RESPONDS: 
GOTO 2500 CART. , CART. PARA. „POLAR, OR POLAR PARA. 7 


INPUT 1,2,3,ОВ 4 


4 PLOTTING WINDOW? 


INPUT XMIN, ХМАХ, YMIN, YMAX 


-1 

1 

-1 

1 PARAMETER/THETA RANGE? 
INPUT TMIN,TMAX | 

-АЖРТ 

ажрт X/THETA EXPRESSION? 
INPUT_X(T) OR THETA(T) 

T/7 Y/R EXPRESSION? 


INPUT Y(X) OR Y(T) OR R(T) 


---------------------------------------------- 


COS I 


(wait about 305) (results plotted) 
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Step 10: Discussion: 


A serpentine, some Lissajous patterns, a Limacon of Pascal, 
and two Rhodoneae are plotted in figures 2-4 and 2-5 . Compare 
your results, keeping in mind the limited graphics resolution of the 
ZX81. Plotting times, which are typically long, have been given 


in each problem. 

If you followed the demonstration carefully, you should have had 
no trouble entering tbe plotting window. However, when you begin 
to experiment beyond the demonstration, remember to enter XMIN,XMAX, 
YMIN,and YMAX in that order (as prompted). Also, remember that 
expressions must contain the appropriate independent variable: X in 
cartesian expressions, T in polar and parametric expressions. 

Values of constants used in expressions may be entered into 
variables used in expressions in immediate mode before executing the 
program. For example, if you have already assigned values to the 
variables YINT and SLOPE, then the expression 


"Y CX) =SLOPE*X+YINT" 


is a valid expression to use. (Remember to GOTO 2500, not RUN, when 
doing this). 


Step 11: Suggestions: 


Using PLOT2 should become easier the more you use it. So try 
some more examples! Use equations and plotting parameters from table 


2-1 and figures 2-4 and 2-5 . Keep track of your plotting parameters 


for each plot. 
There is much room for experimentation in this demonstration. 


Vary the constants in the function of each problem for different 


e‏ حب یت 


results. "Warp" the functions; for example, change a SIN(T) to a 
COS(T), or even EXP(T). 

опе bothersome thing about PLOT2 is that it takes so darned long 
to plot a function. Try experimenting to get faster plots. You 
might try plotting fewer points than 200 (line 2550), or leaving out 
some of the branching or window clipping (lines 2566-2572), or change 
line2574 so that the program doesn't have to think so hard for each 
point. 

If you have ассезз to a higher power machine, try to convert 
this program over just to see how much faster it will run. 

The program шау be split up into 4 different programs, each for 
a different curve representation. Some savings in time may be 
obtained in this way also. 
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Figure 2-4 Plane curves 


Name: Equations of curve: 
a Serpentine Y=X/((X/6)7+1/9) 
b Limacon of Pascal R=2COS(T)+1 , Ө-Т 
с Rhodonea R=COS(4T) , 8=ፐ 
በ Rhodonea R=COS(7T) , в=т 
e Witch of Agnesi Y=2/((X/2)7+1) 
f Eight curve R=2COS(T)* 1+SIN?¢(T) 


Ө-АТМ(ЗІМСТО) 


Т-гапае: 


с-4, 1) 
с-т, Я) 


CTT 


C-*,1) 


Plottin 


(-10,-10), (10,10) 
(=3.5,>2›› С2,29 
(-1,-12, (1,1) 

(8 југа у Ст 12 
(-1,0), (1,2) 


(-2,-2),(2,2) 
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Figure 2-5 Lissajous patterns: 
The equations used for these plots all have the form: 


X(T) =SINCaT+b) 
ҮСТ) =SIN(T) 


where a and b are the constants given in the table below. 
The independent variable, T, was varied over the T-range given in the 
table, and the plotting window in each case was 


(Xmin,Ymin) = (-1,-10 
(Xmax, Ymax)? CTD 


Section: а: р: Т-гапае: 
а . 25 O (-20,201 
b „25 71/16 (-20,20) 
с .25 1/10 (-20,20) 
а .25 1/8 {-20,20) 
e 25 0 (-10,10) 
f anb 1/12 (-10,10) 
g ን 1/6 [-10, 10) 
ከ 5 4/24 (7210,10) 
i #75 о (-20,20) 
j ¿75 1/16 (-20, 20] 
к „#5 1/10 (-20,20) 
1 ‚75 4/8- (-20,20) 
т 1 0 (-4,4) 

n 1 7/8 (-4,4) 

о 1 1/4 (-4,4) 

р 1 1/2 (-4,4) 
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Table 2-1.a Cartesian forms. a,b,c,d = arbitrary constants. 


Мате: (Хх): Х-гапае 
Нурегро1а а/х m ( - 00,00) 
#0 
Witch of Agnesi a/((X/a)**2+1) ( – 20, 00) 
Serpentine Х/( (Х/а) **2+bD**2) ( — өр, 09) 
Quatrix of Hippias X/TAN(axX) (0,4/a) 
Trident of Newton a/X-b«cxX«dxXxx2 ( - 00, со) 
> хо 
Polynomial neo 8j*X**n ( – ор, ад) 
Normal Error Curve a/EXP(* C(axX-b) **2) ( -00, сө) 
сасепагу COSH(axX)/(2xa) ( -00, во) 
Table 2-1.р Cartesian Parametric forms. a,b,c,d - arbitrary constants. 
С = COS(T), S = ЗІМ(Т), T = ТАМСТ? 
мате: X CE ҮСТ): T-range: 
Cross-Curve a/C ህ/= C-*4, T) 
£0,-/2 
Bullet-Nose ажс а/т (-9,7) 
#0 
Lemniscate of Bernoulli ажС/(зжж2+1) ажбжС/ (Зжж2+1) 1-1,7) 
Lemniscate of Gerono axC axS*C (-1,7) 
Lissajous Patterns axs bxSIN(cxT«d) 
Eptrochoid ax(bxC-cxCOSC(bxTO) > ax(bxS-cxSINCbxTO > 
Hypotrochoid ax (b*xC+c*xCOS(b*T)) ax(bxS-cxSIN(bxT)) 
Trochoid ах (T-b*S) аж (1-ржС) ( -со, ор) 
Pedals of a Parabola ахСхх2ж (Тжх2-р) ахбхСх (Tx*2-b) (-q/2, 
1/2) 
Involute of а Circle ах (C+TxS) ах(8-Тхс) ( — од, со) 
Table 2-1.c Polar forms. a,b = arbitrary constants. 
ር = COS(T), S = SINCT), T = ТАМ(ТО 
Name: R(T): Т-гапае 
цітаСоп of Pascal ах (Са р) (-71,1) 
Conchoid of Nichomedes ax (C+b)/C (-9, 1) 
#+1/2 
Карра Сигуе а/Т (0,2ж40 
#4 
Kampyle of Eudoxus a/Cx«x2 L-T.) 
# +ሻ/2 
Folia ахСх (S**2-b) 
Folium of Descartes ахСж$/ (СжжЗ+-9жж3) 
Semi-Cubical Parabola ахТкха/с 
Cochleoid axs/T ( - 00, оа) 
#0 


Table 2-1.4 Polar Parametric forms. a,b = arbitrary constants. 


с = cos(T), S = SIN(T). 
Name: R(T): TOT: T-range 
ت چ ی ی‎ 
Rhodonea axC T/b 

Nephroid of Freeth. аж(8+.5) 2xT (5-1/2;. 
4/22 
Cayley's Sextic axCxx3 эхТ (-Т/2, 
1/2) 
Tschirnhausen's Cubic a/Cxx3 3xT (-9/2, 
1/2) 
Logarithmic Spirals axEXP(CT) T/b ( -00,09) 
Archimedes Spiral axT T ( -09, 00) 
Hyperbolic Spiral a/T T ( -00,00) 

ғо 
Ері Spiral а/с Т/р (-4/2, 
, 1/2) 
Poinsot's Spiral #1 a/COSH(T) T/b ( 200,00) 
Poinsot's Spiral #2 a/SINH(T) T/b ( - (0,00) 

ко 
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Figure 2-6 


The following plots were made on the 7Х81 and TS2068. 

The lower resolution plots were made оп the 2Х81, and the 
higher resolution plots were made on the TS2068 using a 
program similar to PLOT2. 
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2508 сон 
TOUR 

2508 

2504 І, 

EERE 25 

286805 PRINT "СОМТОйВБоЗРАСІНО-СВАА 
МЕТЕВ Т" 

2510 PRINT "INPUT ORDER OF ТО 
2612 INPUT М 

2614 PRINT “PLOTTING МІМОПШЫ?" 
2615 PRINT “INPUT XHIM,XMBEX,vMIM 


1 

EL 

БІ 

53 

TMAX” 
8618 INPUT а 
2620 INPUT В 
2632 INPUT с 
8584 INFUT D 
2626 LET E= 
፳ ZX81, u$ 
2628 LET Е 

R 2:54. us 
2638 LET М 
ЕМ жз For 
дзеів-ну и! 
2632 CLS 
2534 FOR x 
2636 FOR У 
2638 LET Z 
2640 IF ІМ 
(X-B) Е, Ту 
2542 NEXT 
5844 NEXT 
2545 PRINT 
TOURS “15 
2548 PRINT 
СШ 
2550 STOR 


Figure 2-7 Contour plotting program 


PLOTS 


This program plots contours of 3-dimensional surfaces described 
by expressions of the form: 


2 = ЕСХ,У) 


Тһе edges of the contours trace out plane curves іп Х and Y where 2 
has a particular value. The contours are spaced by the use of a 
Contour-Spacing-Parameter, C-S-P. For algebraic expressions, if the 
C-S-P equals the order of the expression (the highest degree in the 
expression), contours should be visible. Spacing may be determined, 
also, by the use of figure 2-8. 


RUN: 
PLOTS: Changes A B C D E F M X Y Z Z$ 
Prompts, іп order of appearance: 
1. "Z EXPRESSION?" 
“INPUT Z(X,Y)" 
2. "CONTOUR-SPACING-PARAMETER?" 
“INPUT ORDER OF Z(X,Y)" 
3. “PLOTTING WINDOW?" 
"INPUT ХМІМ,ХМАХ,ҮМІМ,ҰМАХ" 
ЗТОР: Contours plotted 


Contour spacing and plotting window printed at top of screen. 
Summary of operation: 


Lines 2600 to 2624: Prompt user for expression, contour-spacing- 
parameter, and plotting window. 

Lines 2626 to 2628: Compute step-sizes 

Line 2630: Compute units per contour from С-8-Р 

Lines 2638 to 2640: Plot pixel if it lies on contour band. 

Lines 2646 to 2648: Print units per contour and plotting window. 


47 


т 
=> 
2 
4 
8 10 
er 
~ 
а 
ы 
е 
>= 
~ 
50 
Я Е 
i ч E 
| о E EI 
8 | P 5: 
| D 1 у, LÀ 185 ፻/ 023 2377 65:58 ዘ57' 38 
| ያክ ገ ЧУ ус: ውያ) 
| 4 an: Е 
| 2 
| о 
о 
| в 
| З 
| 21 


| Contour-Spacing-Parameter 


Figure 2-8 Contour spacing vs C-S-P and X-range for program 
PLOT3 
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Demonstration 2.5 PLOT CONTOURS OF SOME SURFACES FROM TABLE 2-2 


Step 1: Enter program: PLOT3 


Step 2 (PROBLEM 1): Plot contours of the elliptic paraboloid with 


equation 
2 = X* + Y? + ХУ 


(-10,-10» 
(10,10) 


Use the plotting window: (Xmin,Ymin) 
(Xmax,Ymax) 


Use a Contour-Spacing-Parameter of 2 since the highest order 


X^ and У? have an exponent of 2: the expression is 2 order. 
gives 20 units per contour, (see figure 2-8). 
Step 3: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
GOTO 2600 Z EXPRESSION? 
INPUT x 
ХжХ+УужУ +Хжу CONTOUR-SPACING-PARAMETER? 
INPUT ORDER ОЕ Z(X,Y) 
2 PLOTTING WINDOW? 
INPUT XMIN,XMAX, УМІМ, YMAX 
-10 
10 
-10 
10 
(wait about 21153} (results plotted). 
Step 4 (PROBLEM 2): Plot contours of the elliptic cone with 
Z = X? + Y2 
Use the plotting window: (Xmin,Ymin) = (-10,-10) 
(Xmax,Ymax) = (10,10) 


Use a Contour-Spacing-Parameter of 1, since the Z expression 
order. This gives 2 units per contour. 


terms, 
This 


equation 


is first 
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Step 5: Run program: 


USER ENTERS: COMPUTER RESPONDS: 


GOTO 2600 Z EXPRESSION? 
INPUT Z(X,Y) 
SQR (XxX-«YxY) CONTOUR-SPACING-PARAMETER? 
____ІМРОТ ORDER OF Z(X,Y) 
1 PLOTTING WINDOW? 
INPUT XMIN,XMAX,YMIN,YMAX 
-10 
10 
-10 
10 
(wait about 6m30s) _ (results plotted) 


Step 6: Discussion: 


The resulting plots have alternate zones of black and white, as 
shown in figure 2-9. The contours, loci of points havimg identical 
Z-values, аге on the edges between these black and white zones. 

Іп step 2, the elliptic paraboloid expression was second order. 
You might think of the variables X and Y as having units, say cm. 
Then, since Z = X^ + Y? + XY, Z must have units ОР см. How about 
expressions like 


Z = aX + DX? + єх? ? 


In this сазе, 2 із third order, taking the constant с to Бе dimension- 
less, i 

Selecting the Contour-Spacing-Parameter to be the order of the 
expression шау or may not give the best plots. To change the spacing, 
use figure 2-8. 

Ав you must have noticed, these plots take tremendous amounts 
of time. (As long as you aren't paying for computer time, who cares?) 
This is because each pixel of the screen must be examined to deter- 
mine its status as either black or white (lines 2634 to 2644). 


Step 7: Suggestions: 


Try varying the Z expressions to give different patterns. To 
visualize the three-dimensional surfaces, label each contour edge, bv 
determining the value of Z for one point and counting up using the 
contour spacing. 

Plot the surfaces using different contour-spacing-parameters. 
In some cases, some interesting "interference" effects occur, 
especially when the C-S-P is small in relation to the order of the 
expression. 

A faster algorithm, but more complicated one, is to start at 
a value of (X,Y) that gives a particular Z-value, and then tracing 
around the contour until the beginning of the contour, or the edge 
of the plotting window is reached. 


Table 2-2 Surfaces defined оп X Y plane. a,b,c - arbitrary constants. 


Name: Z(X,Y): 


Elliptic Cone SQR (ахх) жж2+ (Бжу)жж2) 


Elliptic Paraboloid СажХ) жж2 + (ржу) жж2+сжхжУ 
Hyperbolic Paraboloid Сажх)ужжга-(ржу)жжг 

| SINCOX«iY)| SQR( (COSH(2*Y) -COS(2*X))/2) 
ARG(SINCX+iY)) ATNC(TANH CY) /TANCX) ን 

Plane a*X+bxY+c 


Family of Cassinian Ovals (((X-a)xxg«Yxx2)x(X-«a)x*2-Yxx2))xxb 
Family of. Horse Fetters (ХкжжачУжжа) / (джа) чажУжж2/ (Хжжо+үххоу 


Examples: 
a: b: с: C-S-P: Window: Time: 
1 2. 1 (-10,-10) (10,10) 6m30s 
1 1 1 2 (-10,-10) (10,102 2mi5s 
1 1 2 (-10,-10) (10,10) 2m15s 
.5 (-5,-5) (5,5) 13m - 
Ва 6255,75) (5,5) 12mi5s 
1 1 О 1 (-10,-10) (10,10) 11458 
10 .25 1 (-20,-202 (20,20) 91305 
10 То (-20,-20) (20,20) 31153 
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Figure 2-9 The following contour plots меге made on the 7Х81 
and TS2068. Lower resolution plots were done on the 
2Х81 and higher resolution plots on the TS2068. 
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Chapter 3 SPECIAL FUNCTIONS 


There are many occasions when only complex answers to simple 
questions can be found. For example, here is a simple question: 
What is the circumference of a circle with a diameter of one? 

The solution, Т, is a number which is irrational; that is, no 
ratio of integers will give its value exactly, (although the 
program RATIO, chapter 8, may be used to find an approximate ratio). 

Fortunately for us all, the value of 7 may be determined to 
any desired degree of accuracy by the use of any one of a number 
of techniques. The general idea is to add up a number of simple 
terms to obtain a complex approximation. 

The three functions of this chapter are in the same boat with 
Т : they are all answers to fairly simple questions, but must be 
evaluated in complex ways. 

The first function, erfc(x), is the solution of an "improper" 
integral (one with limits at +/- 00 ): 


co т 
erfc(x) = 2/48 =e at = ца کک‎ 7 at 
- 00 
x ^ 


The integrand is two times the normal error function: 
2 2 
F(t) = 4 exp С 0772/20) 
c / 29 


with m-O and 22.5 


Тһе subroutine ERFCX uses a series expansion to approximate 
values of the complementary error function. Demonstration 3.1 
plots ег#с(х) over O«-X«-2. 

The gamma function is the solution to another improper integral: 


ወ 
Tix = [ሥ‹-‹ሐ 5 x20 


о 
The definition is extended to negative non-integers Бу the recur- 
sion relation which the above formula obeys (see spec-sheet for 
GAMAX). 
When x is ап integer>O, then: 


Го = (Х-1)! = (х-1)(Х-2) ... ።2።1 
where 0! = 1 


The subroutine САМАХ evaluates Гоо using both а series 
approximation and the recursion relation. Demonstration 3.2 plots 
the gamma function and demonstration 3.4 shows how GAMAX may be used 
to give values of factorials. 

Finally, Bessel's functions are solutions to 2፡0 
ential equations of the form: 

x? a? y + х ду + Фа 2 pay = 0 


dx? ах 


order differ- 


Тһе subroutine BESEX finds approximations to the Bessel's 
functions whose values are finite at х-0, (Bessel's functions of 
the first kind), and where Й in the above equation, (Bessel's 
equation), is a non-negative integer. Demonstration 3.3 plots J (х). 
There are many other examples of special functions. The sine 
and cosine functions are probably the most commonly encountered 
ones.  sin(x) is the solution to: 


бу = -y ; yo) = O 


Internal subroutines оп the ZX81 апа TS2068 evaluate sin(x) by 
the use of a series generator. 
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= 
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ЗЕЕ 
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(COMPLEMENTARY ЕН 


[B1 =: 
(51 = 2574. 
(41= 20123 
(Fis SETAE 
23 = 20243 
(11 = BEBTE 
з, дос зе 

v Tel T 

LET xis +i II 


REM AHAS (GRHMA РНЫСТ ION: 
CIM МОЗ) 
LET Wid =: 
LET 467] = 
СЕТ MEER = 
LET 4(5) = 
LET Wisi =, 
LET HECH = 
LET ів) з 
LET 4111 з 
LET iz,ü3 
ШЕТ МЕТНТ 
РОБ Isi Т 
LET ails ! 
MEXT I 

F Жыр 
FoR 121 
LET ЏЕЈ: 
NEST I 
RETURN 
FOR ፲=፳ 
LET Мая 
MEET I 


БЕТ БЕЗЕУ (INTEGER GORDER BE 
FUNCTION ОҒ 15T KIND) 

LET Us-“4K 74 

КЕТ М=1 

FOR Ісі ТОМ 


NEXT I 
LET ШаМ 
РОЯ 151 То ie 
LET мене іде ፲ THR ) 
Е ABS Нноїіє-е THEM GO TO 32 
LET ህ=ህ+ክ 
МЕХТ І 
LET Waele esi TN 
RETURM 


3-1 Special functions subroutines: complementary error 
function, gamma function, and Bessel's functions 


ERFCX 


This subroutine calculates values of the complementary error 
function: 


со. 
erfe(X) = 2A 3 | exec-t* rat ; -сосхссо 
CALL: Х = а 
ЕВЕСХ: Changes I V МО 
RETURN: У = erfc(X) 


Example call: 


LET ERFCX=3000 


LET X=.5 
GOSUB ERFCX 
PRINT "ЕВЕС(";Х; "9 а ";у 
Algorithm (see reference (11): 
-16 
erfc(X) = (s * ах! ) + еггог(х) ¢ |еггог(х) | ‹=3х10 
=1 
a, 0.0705239784 
a5 з 0.0422820123 
аз = 0.0092705272 
а, = 0.0001520143 
as з 0.0002765672 
ав = 0.0000430638 
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ЗАМАХ 


This subroutine calculates values of the gamma function: 


co 
It» з Јоса ; X50 
о 


за мала. = ts 
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Го» - ox $ -00<Х<00 ; XZ non-positive integer 
CALL: X - Argument 
GAMAX : Changes I M V WC) 
RETURN: у = Го 


Example call: 


LET GAMAX=3100 


LET X= 
GOSUB GAMAX 
PRINT "GAMMAC";X;") = "¿WV 


Algorithm (see reference (11): 


1. Let x = fractional part of X 
tm - integral part of X 

2 Toc = і з b,x + error 

= . 
.988205891 
.897056937 
.918206857 
.756704078 
.482199394 
1. 
- 035868343 


NOUS ሠ N 
ዘ 8 ዘ H | 
| | | 

O O G O O O 


ዘ ዘ ዘ 
| 
оо 


pel сосососо 
м 0 
~ 


(бољој коз 


Dae ; 
co» ; 


X+I) 
Іо 


(х) 


2<-Х 
ጊ‹=ጂ‹2 


Х<1 


> 


|еггог(х) | <=3х10 


.577191652 (see Appendix B, Y ) 


-7 


ВЕЗЕХ 


This subroutine calculates values of integer-order Веззе! "5 
functions of the first kind (finite at 0): 


Jp (X) = solution of (Bessel's equation): 
хзағү + X + (X? - Ио“ = O ; Y(X=0) finite 
dx Ж 


where P= integer >-0. 


CALL: М = Bessel's function order СИ іп the above equation) 
X = Argument 

BESEX: Changes I MU V 

RETURN: У з Jy OO 


Example call: 


LET BESEX-3200 


LET N=0 
LET Х-.5 
GOSUB BESEX 
PRINT pU АМА (GM sU ш М 
Algorithm (see references (1), (21, апа (7)): 


1. Тһе series representation of 09 15: 


JyOD =) £-1)! (x/a) 2N 
І%(М-10! 

PET 

2. For a range of values of X, depending on the value of N, the 
series will converge before terms of the séries become on the order 
of the machine round-off error. For J (X), the range of X values 
is about 0<-Х<-15. For larger arguments of X, an asymptotic series 
can be used: See reference (1). 

3. The convergence test used in BESEX is a chegk to see if the 
current term of the series is less than 10 
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Demonstration 3.1 PLOT COMPLEMENTARY ERROR FUNCTION 


Step 1: Enter main program: 


Tid 


i 116) ኸ!” за 


“ከዛ Гаги 


ін 


Әсер 2: Enter subroutines:  ERFCX 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN {wait about 25s) ..(results plotted) 


Step 4: Discussion: 
The normal error function is 


f(x) + 1 ЕХР(-(х-м)?/2 с?) 
c /24 


The function is shown in the figure below for таб апа @?=4 

The complementary error function, erfc(x), represents the Shaded area 
in the figure. If f(x) is the probability distribution function of 
some random variable X, then erfc(x) is the probability that X won't 
take on values inside the range (-x,x). For example, suppose x. 

is the value of a particular measurement of a film thickness minus 

an assumed mean thickness: 


254.” "meas, 7 assumed 


If X has a distribution like that of figure 3-2, then the probability 
that Xi has a value either above x or below -x is diven by: 
РгоБ( | Х| > |х|) + ег?с(х) 


Since а normal error function has area of 1: егрРесО) >= р: 
The demonstration generated а plot ої erfc(x) Рог О‹=х‹=2, and shows 
that erfc(x) slowly decreases from 1 to 0, ав X increases from 0. 


Step 5: Suggestions 
The error function 
erf(x) = 1 - ег?с(х) 


represents the unshaded area in figure 3-2. Modify the main program 
or subroutine to give values of erf(x) for a plot. Also try the 
function 


F(x) = 1 - erfc(x)/2 


which represents the area of the curve in figure 3-2 for -009:<-Х<-х. 
A property of erfc(x) which you may verify by example is: 


erfc(-x) = ~erfc(x) 


Generate а plot for -2‹=Х‹=2. 


Figure 3-2 Normal error function 
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Demonstration 3.2 PLOT GAMMA FUNCTION 


step 1 Enter main program: 
290 RE 

291 LE 

zuz еп! 

293 GÜ 

224 IF 

== IF 

236 PL 

++» ҒО 

#7, Е, 1) 

297 HEX 

Eas вт: 

Step 2 Enter subroutines:  САМАХ 


Step 3: Run program: 


USER ENTERS: . COMPUTER RESPONDS: 
RUN wait about 2m35s (results plotted) 


Step а: Discussion апа suggestions: 


Your plot should be Similar to figure 3-3. The plot limits are 
the same for both plots. The gamma function is discontinuous at non- 
positive integers, as can be seen from figure 3-3. In fact, the 
function is undefined for non-positive integers. 

Try different x-ranges for other plots: the resolution may 
improve for smaller windows. 

For large values of x, Stirling's formula approximates Гю: 


Го» Аз ск/еу“ ኀህ=ግ/፡ (1+1/12x) ; X "lerge" 


Modify the demonstration to compare this approximation with the 
values given by GAMAX. For another variation, modify GAMAX so that 
Stirling's formula is used if X is "large" enough. This will save 
time in the evaluation of | (X), if X is very large. 

Check to see if the subroutine GAMAX obeys the recursion relation 
given in the spec-sheet for GAMAX. “List out a table of values of 
the function for arguments differing by integers. 


Figure 3-3 Gamma function 


-5 
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Demonstration 3.3 PLOT BESSEL'S FUNCTIONS 


3 [በ 
et 


La La a 
£r 


р 1: Enter main program: 


Mee ጠ 


Dae Ricotta = 


ск wala С CA а А 


M CÓ a СИ 


ep 2: Enter subroutines: BESEX 


{Л 


Step 3 (PROBLEM): Plot 2809, іле oth order Bessel's function of the 
first kind, for 0‹=Х‹=15. 


Step а: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT NU=( INTEGER» =0) 
О {wait about 1m15s} (results plotted) 


Step 5: Discussion and suggestions: 


The resulting plot should be similar to that of figure 3-4. 
The oscillatory behavior of the function brings to mind the sine and 
cosine functions. Bessel's functions are like two-dimensional 
analogues to the sine and cosine functions: A round drum head 
vibrates with a radial distortion function that is а Bessel's function 
or a sum of Bessel's functions. By contrast, a string vibrates with 
a distortion function of a sine function or a sum of sine functions. 

The failure of BESEX at large arguments may be handled by an 
asymptotic series. This would increase the complexity of BESEX, but 
could reduce the evaluation times. 

The relation between #6 апа Ji is: 


ЛО = -а J4OD 
ጊ ах р 
In the next chapter, а subroutine to take derivatives із presented. 
Use this subroutine to check out the relationship between Jo and Ј,. 
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Demonstration 3.4 FACTORIALS 


азалы F Ca TO P9 


Ca Ed СА СЯ 


екті 


Step 2: Enter subroutines: GAMAX 
Step 3 (PROBLEM): Calculate 15! = 15*14*13* ... ж2ж1 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN _____________1МРЏОТ INTEGER? =O. 
15 (results printed) 


Step 5: Discussion and suggestions: 
The gamma function is also called the factorial function since: 
n! з Газ) ; · n = integer >-0 


For this example, 15! == 1.3x10!4 

Find the maximum h for which the demonstration program works. 
This may be fairly low. Use Stirling's formula, (see page 62), for 
higher values of n. 

The gamma function may be used to evaluate non-integer order 
Bessel's functions: 


Ју) = rl ЖАУЫ 
£n [Т у+ї+ї› 


Modify ВЕЗЕХ to take advantage of the above formula. 
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Спартег 4 CALCULUS 


In Calculus, there are two operations performed on functions which 
are the inverse of each other. That is, one will "ипдо" the other 
(think of subtraction versus addition). These аге integration and 
differentiation. To perform either one of these operations, the 
numerical approach is to sample the function at a finite number of 
closely-spaced points, and to assume that between each pair of 
points, the function can be well-approximated by an easily inte- 
grated or differentiated function. 

The subroutine DERIV.1 performs a general-order differentiation 
of a sampled function. DERIV.1 assumes that the function is approx- 
imately linear between each pair of points. The first derivative 
of the function, at the point between a pair of sampled points, is 
given approximately by the slope of the line connecting the pair 
of points. 


There is room left at line numbers 4100 to 4200 for a “DERIV.2"- 


which was never written. This subroutine was to "symbolically" 
differentiate functions. That is, instead of chopping the function 
up and approximating the function by a line segment at each point, 
standard rules for differentiation would be used, such as 


п n-1 


а (x) = nx 


ах 


The subroutine INTEG.1 integrates sampled functions. INTEG.1 
uses a general Simpson's rule integration technique. Simpson's 
rule assumes that between pairs of points, the function can be 
approximated by a parabola, which is defined by three successive 
points. INTEG.1 may be used to integrate functions which have been 
sampled at uneven step-sizes. This feature is useful for the inte- 
gration of functions which are rapidly changing in some regions, 
but slowly-changing in others. In addition, INTEG.1 calculates 
error estimates for both even and uneven step-size cases. 

INTEG.2 integrates named functions by the use of а 10-роіпі 
Gaussian quadrature technique. The subroutine samples the function 
ten times, at specific places, and assumes the function is well- 
approximated by a polynomial of degree 19 or less. INTEG.2 is, in 
general, faster than INTEG.1, however, no error estimate is 
calculated by INTEG.2, which may be necessary for many applications. 

The demonstrations in this chapter illustrate the features of 
the three subroutines. For particular applications of the sub- 
routines, some customization may be in order. For example, INTEG.1 
calculates errors based on both even and uneven step-sizes. When 
only even step-sizes are used, some integration time may be saved 
by the elimination of appropriate program lines. 
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Calculus subroutines: Differentiation 


Figure 4-1 


DERIV.1 


This subroutine calculates a general order derivative of a 
sampled function, or a set of data points ((Х(І),У СІ); 1-1,2,3,...). 
The values of X() are assumed to be ordered from smallest to largest 
or largest to smallest, and the derivative at one point, %, із 
calculated for а given value of the index І: 


аб ፕርቋ› 
dx? 
where X = X(I«G) + ха) 


2 


and the sampled function is: Y(I)> = fF(X(I)) , I = 1,2,3,... 
X(1)<X(2)<X(3)<... 
or 
Х(1)›Х(2)›Х(3)›... 


ж 
~. 
ዘ 


CALL: Arguments of function in ascending order, or 
X-coordinates of data points. 
YC) = Function values at corresponding Х(І)'8, or 
y-coordinates of data points. 
а = Order of derivative to be calculated 
I = Index of point where derivative is to be calculated. 


DERIV.1: Changes р) J КУ 
RETURN: v = аб үру ; X = X(I+G) + XO) 
G 2 
Example call: 
LET DERIV=4000 
LET G=2 


DIM X(100) 
DIM Y(100) 


LET I-79 


GOSUB DERIV 
PRINT G;"-ORDER DERIVATIVE OF Y AT "ІХСІ);" = "ју 
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Algorithm: 


1. Define: 
(G) _ а A 
«е а. VOL. ን 
ах 
ጵ. = ХЕ ха 
K,J = 


XCI+K) - X(I+J) 


JX, y 
Дик 3 = Ү(І+К) - Y(I+J) 


2. The general-order derivative is given, in this notation: 


(6) (G-1) _ (6-1) 
зрад G,1 YG-1,0 
x Да 

а,1 ` Ха-1,0 


3. Бог approximately equal stepsizes over Х(І) to X(I«G): 
A А ምቹ 
Васа” Xa-1,0° Аха о 
а 


а. Continuing іп this way: 


у (2 aea! ТИМ ер Да а 


6,0 
АХа,а-і АЖ, ло 


Аха z Па = ረእቾዕ-1., 1 /ቾ6-2.,6 
X 


G, 1 G-1,0 


Aa, o0 


ІМТЕС.1 


This subroutine integrates a sampled function, ог a set of data 
points (‹ጂር፤ንንኛ‹ር፤ንን፲፲=1,2,3,...,ቨዘ3. The values of X() are assumed 
to be ordered from smallest to largest. A generalized Simpson's 
rule summation is used to approximate the integral: 


X(N) 


¥(X) dx 
ха) 


Error estimates аге calculated for both the cases of equal step-size 
and unequal step-size sampling. The trapezoidal rule summation is 
also calculated for estimating errors in "noisy" data integrations. 


CALL: DERIV = 4000 
X() = Arguments of function in ascending order, or 
X-coordinates of data points. 
Function values at corresponding Х(І)'8, or 
y-coordinates of data points. 
N - Number of data samples. (must be. odd) 
= Dimension of XO and ус) 


YO 


ዘ 


INTEG.1: ‘Changes р‹) 53 52 рз ра G I I1 I2 13 J K Q R S T V 


RETURN: 


It 


Simpson's rule summation 
Trapezoidal rule summation 

Equal step-size estimated error 
Unequal step-size estimated error 


5 
T 
Q 
R 


Example call: 


LET INTEG-4200 

LET DERIV-4000 

LET N-21 

DIM X(N) 

DIM Y(N) 

GOSUB INTEG 

PRINT "SIMPSON/S RULE SUM = ":8 | 
PRINT "TRAPEZOIDAL RULE SUM = ";T | 
PRINT "ESTIMATED ERROR = ";Q+R 
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А 


Ақыш ይ ee Раф ~ 


=. ¿n = 


Algorithm (see reference (161): 
1. Define: 
Ateo = X(I+1) - ха» 
Akso = X(1+2) = X(1) 
Ах, у = X(1+2) - X(I+1) 


ፕ'ቾን ۔‎ ፈኞ уска)» 


ax 
2. For the Ма rs Ax, , a Taylor series about X(I) gives: 
X42) У ዎልጾጆእ« 
5(1) = YOOdx = а че ҰО + уа Axe + » 
ха) з 2! 


3. The generalized Simpson's rule is: 


_ (0) (1) (2) л. 3 
SOY = Y А +ኛ Axi Lo T 
2! 


* error(I) 


3 
А (3, Ax _ 

error(1) ү еы к\к. а v I E ን 
Was 


Unequal step-size error 


4. The integral is given bv 


хочу 
Y (X)dX== ува) s I = 1,3,5,...,М-2 
ха) 
5. Тһе estimated error is given Бу 
Error(est) = errori ; l 2 349545552 
1 
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Equal step-size 
error 


INTEG.2 


This program integrates a named function over user specified 
limits. The program uses a 10-point Gaussian quadrature algorithn, 
which is accurate for integrands which may be closely approximated 
by a polynomial of degree 19 or less, over the specified interval. 


RUN: 


INTEG.2: Changes A B I S WO ХО Y$ 
Prompts, in order of appearance: 
1. "INTEGRAND-F(X)" 
2. "INPUT LOWER,UPPER LIMITS" 


STOP: Integral printed 


Algorithm (see reference (1)): 
B 


5 
vooax я») м, | ЕС) + ዞ‹-"1›| 
із! 


F(r) = (B-A)/2 ж 4 (B-A)/2 x r + (B+A)/2] 


eh positive zero of the Legendre polynomial Р, (r) 


> 
2 БЕ стр | 
а - гр TO Уд 


5 
ዘ 
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Demonstration 4.1 PLOT DERIVATIVE 


Step 1 (PROBLEM): Plot the eth derivative of SIN(X). 


Step 2: Enter main program: 


Step 3: Enter subroutines: DERIV.i 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN {wait about 2m205) (results plotted) 


Step 5: Discussion and Suggestions: 
The gra derivative of SIN(XD is "ВІН The resulting plot 
should be of 4 cycles of the sine and its 6 derivative: Х varies 
from 0 to 81, (see lines 324 and 325). That is, the number of cycles is 


п = {Xmax - Хтіп) = 4 
24 


The plot of the derivative is slightly out of phase with the 
original SINCX). This is because the demonstration assumes that the 
calculated derivative is associated with the point (X(I),Y(I)), but 
the subroutine assumes the derivative is associated with the point 


dp,vopo ; where XCI) = xao - ха» 
2 


The difference between Ror апа Х‹1› is small, since the step-size 
is very small. Correct for this shift either by changing the main 
program or changing DERIV.1 itself. Write a more general main 
program to handle any order of derivative, and take care of the 
"shift" problem, which is a function of the order of derivative 


being taken. 
Try different order derivatives of the sine function. To do 
this, change line 328 to 


328 INPUT "ORDER?",G 


Do a number of higher order calculations and get some points for a 
time versus derivative order plot. Use some of the subroutines from | 
chapter 2 to plot these points. i 
The number of cycles for sampling and plotting may be changed by | 
altering program lines 325, 330, and 332. It is interesting to note 
the effect of choosing a sampling such that only a very few number 
of samples are chosen per cycle. Remember that DERIV.1 assumes that 
the function is well-approximated by a piece-wise linear function 
between sampled points. 
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Demonstration 4.2 PLOT FIRST DERIVATIVE 


Step 1: Enter main program: 


Step 2: Enter subroutines: DERIV.1 POINT.2 VINAX SCALE.2 PLOTD 


Step 3 (PROBLEM): Plot the first derivative of cosh(X) for 
-2‹=Х‹=2. 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN INPUT NO. PTS. 

100 INPUT Y EXPRESSION=Y(X) 
(EXP X+EXP -X)/2 INPUT ХМІМ,ХМАХ 

-2 ` 

2 

{wait about 458) . (results plotted) 


Step 5: Discussion and suggestions: 


The derivative of cosh(X) is sinh(X), shcwn in figure 4-2. 


In Appendix А, ("One-liners"5, expressions for these functions аге 
given, in terms of the built-in functions (2Х81 and TS2068): 


Sinh(X) = (EXP X - EXP -X))/2 
cosh(X) = (EXP X + EXP -Х))/2 


совҺһ(Х) is an even function, which is to say that 


cosh(-X) = cosh(X) 


while sinh¢X) is an odd function: 


sinh(-xX) = -sinh(X). 


The first derivative of a function is the slope of the tangents 


to the curve at each point. From figure 4-2, the slopes of the 
tangents to cosh(X) for Х<0 are negative; the slopes of the tangents 
to cosh(X) for Х›0 are positive; and the slope of the tangent to 
cosh(X) at Х-0 is 0. The sinh function indeed shows these properties. 
Some other functions are listed in table 2-1. Differentiate 
Some of these functions by hand and by computer and compare your 
results. 
Modify the demonstration to compile some data points for a 
plot of: difference between calculated derivative and actual 
derivative at each value of X. One might expect the error to be 
greatest for more rapidly-changing regions of the function. Vary 
also: the number of samples, and the function itself. 


Figure 3-2 cosh¢X) and sinh(x) 


TT 
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Demonstration 4.3 SIMPSON'S RULE INTEGRATION (EQUAL STEP-SIZE) 


ИШТЕ 
THE 
с 


CH P. P. ТЕ 


Step 2: Enter subroutines: POINT.2 DERIV.1 INTEG.1 
Step 3 (PROBLEM): Integrate SIN(X) for O< =X< =f 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN NUMBER OF POINTS MUST BE ODD 
| INPUT NO. PTS. 

21 INPUT Y EXPRESSION-Y(X) 

SIN X INPUT XMIN,XMAX 

[8] 

PI (wait about 25s) (results printed) 

---- M29 _______ 5111 1.55. printed) 


Step 5: Discussion and suggestions 
105 result should be close to: 
| занооах = -С08(9) + COS(O) = -(-1) + (+1) = 2 
This 15 the shaded area in figure 4-3. The ratio ОҒ area under the 
Sine function to the area of the block in the figure is the "fill- 


ratio": 


Area under-sine - 
Area under block 


The average value of the function between 0 апа 1 is 


Area under sine = 2 + .64 
Ax 7 


Vary the number of sampled points and note the effect of this 
on the error calculations. Calculate the actual errors: between 
the hand-calculated value of the integral and the numerically inte- 


grated values. Compare the actual error with the estimated error. 


Try some other functions and find out how well. the error estimation 


works for each function. Make a plot of the estimated and actual 


errors versus number of sample points, for each function you inves- 


tigate. Table 4-1 gives some results for the functions: 


EXP(-X), SINOO, 18x? and 8/X*. 


The error calculation used in this demonstration was for equal 


step-sizes. There is also an unequal step-size error calculation, 
used in demonstration 4.4. The unequal step-size errors should 
be small in this demonstration, since equal step-sizes were used. 
Print out this error calculation to see just how small they are. 
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Demonstration 4.4 SIMPSON'S RULE INTEGRATION (UNEQUAL STEP-SIZE) 


Step 1 (PROBLEM): Integrate the function 
У = EXP(-X) 


for Х-0 to the end of 21 steps defined by 


хі = 1.2* AXi 


Step 2: Enter main program: 


Step 3: Enter subroutines: INTEG.1 DERIV.1 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN (wait about 30s} (results printed) 


Step 5: Discussion: 


The unequal step-size sampling is shown in figure 4-4. More 
samples are taken in the region near X-0, where the function is more 
rapidly changing, as compared to the region for larger X. (Plot the 
function's derivative using DERIV.1). 

If enough samples are taken, and from a small enough value of 
X to a large value of X, then the result of INTEG.1 should approximate 
the value of the (improper) integral: 

co 


оо 
Du) «| весах = -ЕХР(-Х) | = 1 
о 


о 


Since the sampling was done with unequal step-sizes, the demonstration 
uses the unequal step-size error calculation. 


Step 6: Suggestions: 


Try this experiment for various functions, including ЕХРС-Х): 
For a range of values of p, sample a function with: 


(.10жр 


ї 


1. Equal step-sizes: Ax, 
(1-1) 


2. Unequal step-sizes: Ах (.1)*p 


To keep the range of integration, R-Xmax-Xmin, relatively constant, 
the number of samples, N, must change: 


1. Equal step-sizes: N INTC1IOR/p) 


E 
2. Unequal step-sizes: Ny = INTCLNC10OR(p-1)+1)/LNCp)) 
Some values of Ng and Ny for R=20 are 
, Ме: Ny? 
1.01 198 110 
1.02 196 81 
1.03 194 65 
1.04 192 56 
1.05 190 4.9 
1.06 188 ዱዱ 
1.07 186 40 
1.08 185 36 
1.09 183 34 
1.10 181 31 


Now, gather some data for a plot of the absolute values of 
estimated and actual errors versus number of samples, N. Do this for 
both the equal and unequal step-size cases. 

This will illustrate the improvement gained by the use of the 
unequal step-size sampling method. 


Figure 4-4 Unequal sampling of EXP(-X) 
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21 INPUT Y ЕХРВЕЗЗІОМ-У СХ) 
ЕХР -Х ІМРОТ ХМІМ,ХМАХ 


Demonstration 4.5 SIMPSON'S RULE INTEGRATION 
CARTIFICIALLY NOISY DATA) 


Step 1 (PROBLEM): Sample the function 

Y(X) = EXP(-X) 
with an equal step-size, for 0‹=Х‹=1. Then add "noise" to the 
sampled data and submit this noisy data to the integration subroutine 
INTEG.1. The "noise" in the example is .5% of the original sampled 
data at maximum. 


Step 2: Enter main program: 


REM CE 
LET ІР 
LET DE 
LET РО 
PRINT чта низ 
20" 
Li 
Er 
E БЧ 


Step 3: Enter subroutines: INTEG.1 DERIV.1 POINT.2 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN NUMBER OF POINTS MUST BE ODD 


INPUT NO. PTS. 


0 
1 


{wait about 25s}(results printed) 


Step 5: Discussion and suggestions: 


This demonstration presents an error estimation for integration 
of data corrupted by fairly large values of "noise". The "noise" in 
a real-world situation could be introduced by measurement inaccuracy, 
or fluctuation in the measured quantity. The error estimate 
in this demonstration is not a very good estimate of the actual error 
but at least gives a "ball-park" number. 

Try larger values of "noise" and different functions. Also 
try different "noise" calculations, by modifying line 397. 


Table 4-1 Some results from integrations using INTEG.1. The 
functions were sampled N times over the interval (A,B). 
H is the step-size. Estimated error and actual errors 
are given as well as the ratio of the estimated to actual 


error. 
Time to Error Error 
Integrate: (Actual): (Estimated): Ratio: 
3s 1.1 „14 · 125 
125 „040 . 0098 . 242 
235 . 0058 . 0023 -401 
728 .000071 . 000052 .723 
SINCX) 7 0 Т 0.52 55 . 00086 - 0011 1.26 
21 0.16 195 .0000068 .0000070 1.00 
35 0.09 335 .00000081 .00000082 1.00 
105 0.03 1003 .200000009. 000000009 1.00 
18፪7 9 0 1 0.12 7S .0015 .0016 1.12 
27 0.04 258 - 000013 „000014 1.06 
вих? 11 1 2 0.10 98 . 0001 . 000065 .65 
зз 0.03 315 .00000098 .00000084 .86 
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Demonstration 4.6 


Step 1: Enter program: 


Step 2 (PROBLEM): 


F(X) 


= X + X? + X + 


for О<=Х<=10. 


Step 3: Run program: 


GOTO 4300 


X*X*X+K*XK+K+1 


0 
10 


USER ENTERS: 


GAUSSIAN INTEGRATION (POLYNOMIAL) 


INTEG.2 


Integrate the function 


1 


COMPUTER RESPONDS: 


INPUT INTEGRAND=F(X) 


INPUT LOWER, UPPER LIMITS 


(results printed) 


Step 4: Discussion and suggestions: 


The function 


F(X) 


= x? + х2 + Х + 


1 = (X + 10(Х% + ጊን 


is а polynomial of degree less than 19, so ІМТЕС.2 gives ап accurate 
result for its integration, which 15: 


ю 
| rax " (е pK 
vo 


4 3 


2 dd = 
& y x) | - 2893.33 
2 e 


Integrate this function over different limits, and check the 
accuracy for different limits. Try higher degree polynomials. 


Use PLOT2, 


from chapter 2, 


its integration. 


for a clearer picture of the function and 


Demonstration 4.7 GAUSSIAN INTEGRATION (NON-POLYNOMIAL) 


Step 1: Enter program: INTEG.2 

Step 2 (PROBLEM): See how well INTEG.2 integrates the function 
YCK) = (X-1)/LNC(X) 

for О<=Х<=1. 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 4300 INPUT _ INTEGRAND=F(X) 
(X-1)/LN X INPUT LOWER,UPPER LIMITS 
0 

1 (results printed) 


Step 4: Discussion and suggestions: 


Even though the function is not a polynomial, the result is 
close to the exact result: 


' 
Је (X-1)/LN(X) = „75 
о 


Apparently, the function is well approximated by а polynomial 
of degree 19 ог less, at least in the interval 0<-Х<-1. 

Try different limits and compare the results with hand-calculated 
values. Try integrating some of the functions in table 2-1. 

Reference (1) gives different versions of the Gaussian quadrature 
technique. INTEG.2 may be improved if a higher version is used. 
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Chapter 5 MATRIX OPERATIONS 


To see the forest despite the trees is a virtue well documented. 


For example, if we are given a set of equations to solve: 


> = 
Зх, + 2х, 5 
2х. + x, + Хо = 7 
x, + Ха 2 


it is possible to think of the problem опа higher level: the 
transformation of a vector by an operator matrix: 


i 5 
& = x and ë = 7 
2 2 
Хз ее 


Matrices and vectors obey general rules, so that without look- 
ing too deeply into the details of the problem at hand, certain 
conclusions may be drawn. For example, if A's determinant is 0, 
the problem may have no unique solutions. 

This higher-level approach to the problem is an example of 
"chunking", and its advantage in this particular case is that it 
gives a method whereby an infinite number of similar problems may 


be solved: 


1. Call a subroutine to invert А. 2 
2. Call а subroutine to multiply С by A's inverse. 


The two subroutines in this chapter accomplish these two steps. 
MXINV inverts a matrix, апа MXMUT multiplies matrices. Іп both 
subroutines, the matrices to be multiplied or inverted are named 
in a string before calling the subroutine. This is the same 
approach as the vector operations subroutines. 


Several bench-mark inversion times are given in FIG- 5-2.. 
These are relatively long, and the reader may be able to find more 
efficient algorithms for matrix-inversion. It would also be a 


challenge to write the subroutine in assembly language, for 
much higher speeds. 


28088 БЕМ HAIN (HATRIX IMMERSE |p 
ETERHIHRHT! 

seas DIM LUiNH,2zhN! 

5004 DIM МІМ,М) 

5005 FOR Іші To ከ 

2008 FOR 1-1 TO N 

5010 LET |[|(፲.,ህ) svAL xs 

2012 LET Url, J+H) =8 


5014 NEXT ш 
5016 LET ШІІ,ІзМісі 
2818 НЕХТ I 
5020 LET Us 

5022 РОВ Iz 

5824 LET ES 
Бәсе FOR J= 


5042 50 T 
5044 FOR 
5846 LET 
5045 LET 
5058 LET 
5053 LET 
5854 МЕХТ 
SESE РОВ 


sas: 
55652 IF J 58 
58 

S062 FOR 

5064 LET 

5066 NEXT 

505 МЕХТ J 

5070 МЕХТ І 

5878 БОБ Ізі TO H 

5074 FOR 1-1 ТО м 

5076 LET ПІТ, Чі ат, фам) 

5873 NEXT J 

5280 NEXT І 

2822 RETURA 

5108 REM НХМОТ (HBTRIX-HMULTIFLIC 
ATION) 

5182 DIM имал) 

5104 FOR Іші То H 

5186 FOR К=1 TO L 

5128 FOR (ші ТӘНІ 

5118 LET ШІІ,кізімі,кІ жоны ка 
5112 NEXT J 

5114 T É 

5116 МЕХТ Т 

5118 RETURN 


Figure 5-1 Matrix operations subroutines: 


and matrix multiplication 


Matrix inversion 
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МХІМУ 


This subroutine finds the inverse and determinant of а matrix 
named by the string X$. If the matrix is singular, (determinant = 0), 
an error message is printed. 


CALL: N = Dimension of named matrix 

X$ = Name of matrix, indexed by I and J 
MXINV: Changes I J K L М U 00,» VC,D 
RETURN: У(,) = Inverse of named matrix 

U - Determinant of named matrix 


Example call: 
LET MXINV=5000 


LET N=3 
DIM АСМ, ND 


LET X$="A(I, J)" 
GOSUB MXINV 
PRINT "DETERMINANT = ";U 
Example of algorithm (Gaussian reduction, see reference (1110: 


1. Invert the 2x2 matrix А: 


~ 
А = 1 4 
з 2 


2. Lines 5000-5020: Load the left-hand side of the partitioned 
matrix ሾ with X and the 885 of Ü with the 2x2 identity: 


- і 
) з 1 4 і 10 
з 2% 0 1 
' 


Initialize the determinant variable: 


Ú = 1 
~ 
з. Lines 5024-5034: Find largest element іп column 1 of U, and its 
index, K: 
K = 2 
L = UCK,1) = 3 


д. Lines 5036-5042: Update determinant variable: 


U = Ux(-L) + -3 


10. 


ll. 


12. 


Lines 5044-5054: Exchange rows K and 1, then divide each 
element of the new first row by L: 


ы Й 
Ш-|і 2/3 ! 0 1/3 
1 4. 1 1 О 


Lines 5056-5068: Subtract row 1 elements from row 2 elements, 
column by column: 


0 =|1 о/з ' O 1/3 
0 10/3 ‘ 1 -1/3 


Lines 5024-5034: Find largest element in column 2 of ችኮ. апа 115 
іпдех, К: 


К з 


L (K,2) = 10/3 


Lines 5036-5042: Update determinant variable: 
U = Ож = -10 


Lines 5044-5054: Divide each element of the second row by L: 


О 
1 
0 1 | 
1 


Ú «(1 2/3 
3/10 -1/10 


0 1/3 | 


Lines 5056-5068: Subtract 2/3 times row 2 elements from row 1 
elements, column by column: 


і 
бо«Ї1 o ! -2/10 4/10 
0 1 ! 3/10 -1/10 


Lines 5072-5080: Load V(,) with the RHS of Й: 


ፍታ 
У = |-.= „А 
«З "dl 


Тһе result is 
XV 
det(X) = ሀ 


89 


—— Ia... a —— ....- 
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MXMUT 


This subroutine performs matrix multiplication between two 
matrices or vectors named in the string variable X$. The result 


is put 


CALL: 


MXMUT : 


RETURN: 


into the matrix 


Х$ = String 


ሀ(.›. 


containing the names of two arrays or vectors 


to be multiplied: 
1. Separated by а "x" 


2. The 

Ij 
3. The 

J,K 
Number 
Number 
Number 
L = Number 


N 
M 


шин 


Changes I J 


left-hand-side matrix or vector is indexed by 
or J ` 
right-hand-side matrix or vector is indexed 
OF J 

of rows in left-hand-side array 

of columns in left-hand-side array 

of rows in right-hand-side array 

of columns in right-hand-side array 


к ሀ(.› 


UC,) = result of matrix multiplication 
0 has dimensions U(N,L) 


Example call: 


LET MXMUT-5100 
LET N-4 

LET Ме5 

LET L-6 

DIM АСМ,М) 


-DI 


M B(M,L) 


LET X$="ACI,J)*BCJ,K)" 
GOSUB MXMUT 


Example of algorithm: 


1. Calculate the 


зі 


Тһе. 
Тһе 
The 
The 


The 


product of matrices x and É: 


-ХжЖ- (1 чук | 5 "| 


Іп 1115 сазе: > 
number of rows іп А: N 
number of columns in A: М 
питрег of rows іп В: 

number of columns in W: L 


4 6 8 


1 
2 
3 


result matrix will have dimension U(1,3) 


3. Each element оғ Ü is 


of X with the column 
П(1,12 = «(1 2ን›, (3 
061,22 = «(1 2), (5 
UCL, 89 = COL 2),(7 


4. The result is 


ሸ 


(11 17 23) 


the Euclidean dot-product of the row vector 


vector of Ж, 


42) 
6 
а)» 


1х3 + 2x4 
15 + 2x6 
1x7 + 2x8 


(see 


demonstration 1.2): 


11 
17 
23 
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Demonstration 5.1 MATRIX INVERTER 


Step 1: Enter main program: 


418 REM СЕМОБ,1 
411 LET ዘዛጻ ፲8ህ=2ቿፎቿፎ 
412 | ሮ 
412 
414. 
415 їні 
4185 1 
ф17 Tz 
415 zm = 
CROLL ; 
419 PRINT "аі"; dsl; то" 
420 INPUT AIL, Ji 
221 PRINT AIT, Ji 
422 MEXT J 
423 REM #++ FOR USE зі = 
CROLL 
424 PRINT ' ' 
425 MENT I 
ШЕ CLS 
ЧЕТ LET 
425 за 
чая PRI 
422 FOR 
431 FOR 
432 PRIM 
š ace 
ВЕТ 
ТЕ 
438 
PRINT 


ба ጀማ 


RRRERPRRORAPREC 
= ЕН 


PMI GE C I rara 
BEMIS TMN Ра 
бан 


їй 


Step 2: Enter subroutines: MXINV 


Step 3 (PROBLEM 1): Find the inverse of the matrix 


ያ 
А = |1 2 
0 1 


Step а: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
UN INPUT DIM 
2 INPUT MX ELEMENTS: 
———— F С IR? _ ашады. 
+ А(1,2) =? 
; 2 A(2,1)=? 
| 0 A(2,2)=? 
1 (inverse, determinant printed) 


Step 5 (PROBLEM 2): Find the inverse of the matrix 


Xe 


OOOO 
сооню 
O OPF NO 
Or NOT 
мора 


Step 6: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN | 

5 INPUT MX ELEMENTS: 
AC1,1)=? 

1 А(1,2) +? 


enter matrix elements row by row 


"———— E Віз በሼ 
1 (wait about 17s) (results printed) 
PRESS "ENTER" FOR MORE 
(enter) 
(wait about 4s) (results printed) 


Step 7: Discussion and suggestions: 


The matrices in step 3 and step 5 were both invertible: neither 
one of their respective determinants was zero. The two matrices were 
constructed so that the diagonal elements (ል(1,1), „АСМ,М)) were 
all 1 апа the elements below the diagonal were all “0. This guaran- 
tees that the matrices can be inverted, since the determinant in 
these cases is the product of the diagonal elements, or 1. Similarly 
constructed matrices, with various dimensions were inverted, and the 
inversion times are plotted in figure 5-2. The time required for 
the inversion of a matrix follows the relationship 


toc 2-07 


where М is the dimension of the matrix. 

Check the results of this demonstration by multiplying the 
original matrices by their respective inverted matrices. Use the 
matrix-multiplication subroutine MXMUT. The multiplication should 
give the identity matrix of the same dimension as X: 


І з 1 O^ wes O 0 
O 1 0 0 
0 0 1.0 
0 0 O 1 
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Find the inverted matrices of some: 


1. Symmetric matrices: А(Т, 7) = АСІ,І) 
2. Identity matrices: ACI,J) = O if 100] 
1. qf ЕЈ 
з. Diagonal matrices: ACI,J) = O if 1<>J 

1608 | і ИШТ ТТТ 

во Ë ЕН НЕ: . 

d: НА 


20 


a 


| ТЕ | | 


тов | 
8 ШШШ ШИН 
7 . 
4 || ЇН 
m ШІ ШІП 
u 5 E 1. ІШ ШЕН ШШ 
= 4 |. ЧУДИ | a | | 
2 m re v Ана ревно ІП 
" H: || АШ | / 11 | 2. | 
__| | | | | | | | || | 
25 P | | | ^ | a А ЕШ 
Eoo 
ШШ HHH | | 
— 1.1 
' 1 2 3 4 5 6 78910 
DIM 
Figure 5-2 Inversion times vs dimension for matrices of same 


ва form аз those of demonstration 5.1 


Demonstration 5.2 MATRIX-MATRIX MULTIPLICATION 


Step 1 (PROBLEM): Multiply matrix А by matrix É: 


— 0.6 1.0 1.4 1.8 2.2 
ВА = 0.8 1.2 1.6 2.0 2.4 
1.0 1.4 1.8 2.2 2.6 


fO N م س‎ 
онов e م‎ 


Коби ди Анти 


፡የበ[ባቻ ШІЛ ететі) 
x 

እ 1 

A 


TOT ECT ET AAG 


ЕЗІ 


Step 3: Enter subroutines:  MXMUT 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN (results printed) 
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i Demonstration 5.3 MATRIX-VECTOR MULTIPLICATION 


Step 1 (PROBLEM): Multiply vector X ру matrix В: 


оғы 
OMOFfN 


CH en 
nimimm 


пог ERE DDE т 
~ 
፳=2-144-42 


FOR 


РРЕРРЕРРРРЕРРРЕРРРЕРРЕР 
i J TF Ca ከ” 53 un C — J hi ERE С TO e ҒӘ 
т 
e 
n 


MC i IC un wn UPC VOR O С ж OG CO CO PO Cu Cr Т 


Сі 
їй 
+ 
ር; 
ዝ 


Step 3: Enter subroutines: MXMUT 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN (results printed) 
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Demonstration 5.4 TRANSPOSED УЕСТОВ-МАТВІХ MULTIPLICATION 


Step 1 (PROBLEM): Multiply the matrix B Бу the transposed vector x" 


፻ ጃጀ (1 1 1) 0.6 1.0 1.4 1.8 2.2 

0.8 1.2 1.6 2.0 2.4 
1.0 1.4 1.8 2.2 2.6 

Step 2 Enter main program: 

ша REH 

Sal LET 

SAZ LET 

283 LET 

284 LET 

505 DIH 

SUS DIM 

507 FOR 

508 FOR 

Е 

5: 

Step 3: Enter subroutines:  MXMUT 

Step 4: Run program: 

USER ENTERS: COMPUTER RESPONDS: 

RUN (results printed) 

Ue ue ts printed) | 


Step 5: Discussion: 


Demonstrations 5.2, 5.3, and 5.4 of the subroutine MXMUT show 
how any arbitrarily dimensioned matrix may be multiplied by another, 
as long as the number of columns of the left-hand- Side matrix is 
equal to the number of rows of the right-hand-side matrix. The 
two subroutines, MXMUT and MXINV are used in the following two 
two chapters to solve a set of linear equations: 


ARSC 


Y 
where X is the solution vector, А is a square coefficient-matrix, апа 
с is a constant vector. The solution is 


97 


The implied multiplication is of а vector by а matrix, and ДРА is 


the inverse of 
For example, the following set of equations: 


xT 


can be written in matrix form: 


TIE 


The solution is 


Step 6: Suggestions: 


Demonstrations 5.2, 5.3, and 5.4 can be made more general by 
the addition of a matrix-element епігу portion. This might be a 
separate subroutine, MXELM, which could also be called by any other 
program or subroutine to enter elements into arrays. Once you 
bave this working, try multiplving several vectórs, X, by the matrix X: 


Нов [Е 


ТЕ you find a vector X which when multiplied by X gives a multiple 
of itself: 


БЕШЕ 


then you have found ап "eigen-vector" of the matrix X, and K is its 
corresponding "eigen-value." 


Chapter 6 ROOTS/REGRESSION 


This chapter contains three programs and subroutines to solve 
two common numerical analysis problems. The first problem is to 
find a solution of a set of equations which are not necessarily 
linear. The second is to find an equation to "fit" a set of data ፲ 
points. The three programs and subroutines, ROOTN, REGRP, and РЕСЕМ 
аге included together іп this chapter, because their structures 
are similar to one another, and they all make intrinsic use of the 
subroutines MXINV and MXMUT, from chapter 5. 

The program ROOTN solves a general set of equations: 


f (Ха, „ха - 0 
f 
NO, .. ХИ) + 0 
The user must first supply ап initial estimate of the solution 
vector (x .. Xu), and the program iterates until а (not necessarily] 


unique solution» is found. 

The subroutines REGRP and REGRN find a set of parameters 
(w,) for an equation to fit a set of data points: ((X,,...,X,,V2). 
Thd resulting choice of parameters minimizes the error between the 
fitting-function and the data. In both subroutines, the form of 
the fitting equations must be known beforehand: 


In REGRP: የ(ጃ,ወ)ን = Wi FIOR uu X +... + м, f. СХ pees Xp) 
where the basis functions ፻..- -.'ከ are simply functions of 
Хр]... „Хе, зисп ав exp(x. жогу. 

Іп REGRN: Е(Х,М) = any General function of x,,...,X. апа 


1 Е 


parameters (м.). 

REGRP calculates the parameters (w,) without any iteration, 
while REGRN requires both a user-suppliéd initial guess for the 
parameters, as well as several iterations to give the solution. 

The reader and user of ROOTN, REGRP, and REGRN should have no 
trouble in finding real-world problems which can be solved by one 
or more of these programs and subroutines. This is their best test, 
and many improvements to the routines may be found in this way. 


| 
| 
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= 
ті 
сі 
is 
a 


[НЕНТОМ E 


Sana REM РБОВЕАА! 


де 
Lee 


EU 
NO, EQNS, VARS, 


і Beas INPUT 
ваја СІМ 
5012 DIM 
6014 ‘CIH 
6016 DII! 
ваја СІМ 
вага ЕО 
ВАРЕЗ PR 
ИКАТ 
524 IH 
5028 FR 
бога гы 
ваза CL: 
ваза МЕХ 
8034 LE 
БОЗЕ PR 
ваза PR 


Е 
P 
F 
d.n 
База FOR 
5242 LET 
52844 LET 
geiris,ug8 
50465 PRI 
Gl 
Б245 MEST J 
2850 PRINT “----------- 
баба PRINT "CONTINUE? 
BGS4 INPUT хо 
555 IF X$-2"M" THEN STOF 
вага с! + 
Есей Far 
ЕВЕ For 
EMEL LET 
2066 LET 
SG@68 LET Ғфіші-ҒБіші) 
"DO URL 
ЕР LET хі 
Бате HEST Е 
2074 NEST J 
S076 LET ка АТ. gi" 
вата GO SUE MINU 
вага LET ዛ=ክ 
Saa LET | зі 
8084 LET x=" VJ, N Ftd" 
555 GO SUB HERT 
BOSS ЕП зі і 
База LE і 2) 
5092 HE: | 
2894. LE +1 
ease со BIE 


Figure 6-1 Root finding program 
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ROOTN 


This program finds a solution vector, not necessarily unique, 
to a system of equations 


Е С Хр Ху) = 0 ; 

Fy GG »ጆ>›---» Xy? = 0 E 
The program uses a general Newton's қ The user supplies ап 
initial guess for the solution vector, Ху), and the program 4 
iterates to a solution. The number of ieee is equal to the 
number if independent variables, (X.). At each iteration, the 


current value of the solution vector is printed. A root is found 
when the set of function values evaluated at the root is zero. 


RUN: 


ROOTN: INITIALIZATION: 
Prompts, in order of appearance: 
1. "INPUT NO. EQNS./VARS." 

2. "INPUT FI(X(150,...,X(ND2" 

3. "INPUT INITIAL X(I)" 


(where I = 1,...,N) 


PRINT OUT X(1), 


PROMPT: 
"CONTINUE?" 


NEWTON ALGORITHM 


»X(N) апа F1,...,FN 


STOP: 
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Algorithm (General Newton): 


Given a current grs iteration estimate for the solution 


Ñ. = (K. нау 3 
I 11 Ny 
| of the set of equations 
4 
3 =. w 
| F CK) = 0 
Py CX) з 0 


the пеш estimate, x , is 
I+1 
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дя) 


0% 


де“) 


0% 


Demonstration 6.1 SYSTEM OF LINEAR EQUATIONS 


Step 1: Enter program and subroutines: ROOTN MXINV MXMUT 


Step 2 (PROBLEM): Find the solution vector to the following system 
of linear equations 


1 2 з а Х(1) зо 
o i 2 3S X(2) = | 16 
0 0 1 2 X(3) 7 
0 0 0 1 X(4) 2 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 6000 INPUT NO. EQNS. /VARS. 

n INPUT F1(X(1),...,X(N)) 

Х(а) +2жХ(2) +ЗжХ(3) + АжХ(А)-30 INPUT INITIAL Х(1) 

1 INPUT F2(X(1),...,X(N)) 

X(2)+2*XK(3)+3*K(4)-16 INPUT INITIAL X(2) 

1 INPUT F3(X(15,...,XCNDD 

ዚ(3)+2።#ጂ‹ዱ›-7 INPUT INITIAL X(3) 

1 INPUT F4(XC1),...,% (ND) 

Х(42-2 INPUT INITIAL Х(4) 

1 (iteration О printed) 
CONTINUE? 

Center} (iteration 1 printed) 
CONTINUE? ч 

(enter) (iteration 2 printed) 
CONTINUE? 


N (program stops?) 


Step 3: Discussion and suggestions: 


The system of linear equations problem is easily solved by the 
Newton algorithm in one or two iterations. This is because the 
algorithm essentially linearizes the set of equations at each 
iteration. For example, when there is only one independent variable, 
a general problem might look like that of figure 6-2, where the root 
of F(X) is to be found. In the figure, the current estimate is 


X.. The program finds the next iteration estimate, by: 


I Яра 


1. Find the tangent to F(X) at X.. 
2. Krad = X-intercept of the tangent line. 


Find the number of iterations needed for various sizes of the 
coefficient matrix. Try some matrices that have elements with very 
large and very small values, to see if more than 2 iterations are 
required to find the solution. 

Draw a picture of the general 2-dimensional root finding 
problem and compare it with the 1-d case shown in figure 6-2. 
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Figure 6-2 Root finding problem for 


-10 
-10 0 10 
X1 


Figure 6-3 Root of demonstration 6.2 


TANGENT 


one 


independent variable 


Demonstration 6.2 SYSTEM OF NON-LINEAR EQUATIONS | 


Step 1: Enter program and subroutines: ROOTN МХІМУ MXMUT 


Step 2 (PROBLEM): Find a solution of the following system of 
non-linear equations: 


EXP(X,) = 1 
ү 
ቾ1 + X, = 2 
with initial guess: Xi - Ха = 1 
Step 3: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
GOTO 6000 INPUT NO. EQNS./VARS. 


— –_____МРЏТ NO. EQNS./VARS.  — 
2 INPUT Fi(X(1),...,K(N)) 
DAE X(2)/X(C125-1 INPUT INITIAL X(1) 

INPUT F2(X(1), „Х(М)) 


xG 2+X(2 INPUT INITIAL X» 


(iteration O printed) 


CONTINUE? 
Center) ^ iteration 1 printed) ^ 
CONTINUE? 
keep iterating! 
Center) (iteration 7 printed) 
CONTINUE? 
N (program stops) 


Step 4: Discussion and suggestions: 


In this problem, the common zero of the two surfaces F1 and F2 


is found. From F2 - Xi * Ха - 2, the zero must lie оп the line 


Ха = 2 - X, 
From F1 = ЕХР«Х.2/Х, - 1, the zero must also lie оп the curve: 
Хо = LN(X) 


The common zero is the intersection of these two curves, shown in 
figure 6-3. 

Try different initial guesses, some from each quadrant of the 
Х,-Х, coordinate system, shown in figure 6-3. Find out which values 
converge most rapidly to a solution, and which, if any, never converge; 
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Demonstration 6.3 NON-LINEAR EQUATION 


Step 1: Enter program and subroutines:  ROOTN MXINV MXMUT 


Step 2 (PROBLEM): Find a solution of the following non-linear 
equation 


COS(X) = Х 

with initial guess: X = 1 

Step 3: Run program: 

USER ENTERS: COMPUTER RESPONDS: 

GOTO 6000 INPUT NO. EQNS./VARS. 

1 INPUT Е1(Х<120,...,Х(МОО 

COS X(1)-XK(1) INPUT INITIAL X(1) 

1 (iteration О printed) 

.. CONTINUE? 

{enter} (iteration 1 printed) 

CONTINUE? 
keep iterating! 

{enter} (iteration 4 printed) 
CONTINUE? 

N (program stops) 


Step 4: Discussion and suggestions: 


The one-dimensional problem in step 2 is shown in figure 6-4. 
From the figure, there are many "local minima” in the function 
COS(X) - X. That is, the value of the function at a local minimum 
is less than the values of tne function in the immediate vicinity. 
This may cause trouble for ROOTN, since the algorithm requires 
the X-intercept of a tangent line at X, be closer to the root than 
х However, аз long as the initial guess is within some range 
close to the root, the requirement is met in this problem. 
Determine the range of values of initial guess which give a 
successful convergence of the program. 

There are also problems where more than one zero exists, In 
this case, the initial guess is also critical to the outcome of the 
program. Try some problems with many zeros, such as F(X) = COS(X). 
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ROOT 


cos (x)- x 
o 


-10 


Figure 6-4 Root of COS(X)-x 


LOCAL MINIMUM 
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БІЙ БЕН БЕШЕР (PLS REGRESSION) 
віза DIM ZiM, NI 
2184 СІНО МІНІ 
=5ኋፎ5 DIM СІМ) 
ва FOR Іші To Е 
ла FOR dei TO አ 
LET Grid) + НА. Fidi 
HEAT М 
FOR Ј=1 TG ከ 
FOR Қа4ғ1 To H 
LET ZiíiJ Rl =Z id КЪВ.) етік) 
LET ЕК, И а dl, КІ 
МЕХТ К 
LET ran Е, زل 5+ اا‎ Ги: 
LET ПІІ si EE إل + !ا‎ g£ r 
HEAT < 
ЧЕХТ І 


LET Хайді, 1" 
GD SUB мхтнъ 
LET Мам 


LET Ea ITT LN ғ 
GO SUB М>МЫТ 

FOR Ја] ТО M 

LET iid) lita, ኋን 

МЕХТ 4 

RETURN 


Epp EP ES PS IP ES PS S ES EE ES E P E ра рз EE 


En Cr OS En CI Ch Ch ER ER стис CR Ef C с ски ЄВ Cf сл Єр Gf E ги 
OPP ዮ ዮዮ ርነ CATA CD O ክነ ከ፲ NO NÛ DI pz EP ES I 


ከ) Gy Cp (n p= ከ)6 0) እዮ [O EJ Oj Chi P TO C3 бй ጦዮ Па 


Figure 6-5 Linear least-squares regression subroutine 
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REGRP 


This subroutine performs a linear least-squares regression of 
a set of data points 


COX езеро ГО a er 2ኛ -3ታ 
11 Е, 1 1р Ep P 
to a function which is the linear sum of "basis functions", (ኛ]) 
BOR) = ED ፈጀ t ሣህ P 
where ¥ = Ks. cook) 


For example, а suitable basis function is Е ርጂ. „ХО = 1 + Xi * 2X,, 


The weights w .,W, are found, such that tnê efror 


= ҺЕ 2 
ЕГГОГ = (FOX) YQ) 


is minimized. 


- Number of data points 


CALL: Р 
Е = Number of independent variables: Х pee ХЕ 
М = Number of basis functions: ЖОНЛЕ ғ 
= Number of weights: w pene We : Р 
X(,) = Array of independent variable values of data points | 
YQ) = Array of dependent variable values of data points 
REGRP: Changes GQ) I J K LM ሀ UC) У(,» WO X$ Z(,) 
RETURN: WO) = Array of weights for each basis function 
Algorithm (linear least-squares regression): 
Given P data points 
(<%,,Ү1),..., OV YU) $ X, = сх, е ХЕ 2 


І І 
апа М basis functions, (Е OD); the parameters, tw), for the 
regression function with tne least error 


F(X) = Е, OD (X 
= Wy 1 + ele + Wy Fy X) 
are given by 
-1 
wi «Е „Е? а: SF Fy? CY, Pup 
Wy Pye Fp? aus CF CEN SY FA 
= ә = E F X > 
where (Е БЕК? кұқық) 3 <Y,F > Ж jt I 
те Та, 
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Н Demonstration 6.4 GENERAL POLYNOMIAL, REGRESSION 


| Step 1: Enter main program: 


| 530 REM РЕМОБ,4. 
| 531 LET БЕСБРІВ 
| 538 LET HXINU-EBOD 
| 533 LET 11ጻ[ህፐ=5 
| 5: NT “INPUT NO. D 
| 535 INPUT P т Mere С 


INPUT М 
PRINT “INPUT мо. 


== 


== 
JAR 
за 
Фа DIM FR 
41 | 
42 ( 
43 = 
44 


¿s= 
4 
= 


ІЛ Cn ЕЛ (በ ЕЙ 


AR CEA ЕЛ Ея Ссс. (በር 


ХІМ, 1 1” 
INPUT Fir) 
፲ 


МЕХТ 


ሥ 
„> 
«о 


нт poc 


መም 
a 


г 
ОН 
T 


ባበዬቶዬዮዬ 


Ма а 688) ER ~ ТИЈ 


лат 


Vn cn On on rr 


E E — шо-итил E Ca ff 


mean пелени 
4 


RR Dine 


፲ሀ8(8 =፤በ፤በሂቫሂበ 
ወ 
mmc mi 


r T 
prn 


ела = ሂካ on on 


=4፪ቨ ፪ሽ cn en en — c enm 


Tre ru c ~ С n F eC 
т 
ا‎ 
፲) 


OR 121 ТО N 
LET M=HLUAL Бари suru 
МЕЗТ .! 


= 
5 
= 
= 
= 


“4-4 


Step 2: Enter subroutines:  REGRP МХІМУ MXMUT 


‘Step 3 (PROBLEM 1): Given the data in the table below, find the 
coefficients (W(I)) of the least-squares fit to the function 


F(X,.X,) = МО) + W(2)*EXP(-X,) + W(3) *EXP(-X,) + W(A)*EXP(-X, ።ጁ,) 
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ҥч 
ж 
= 
ም 
еч 
Бед 
ж 
м. 
м 
— 
E 
~ 
= 
хм 


1 1 1.00 6.2 
2 + 1.25 6.0 
3 4. 1.50 le 8 
4. 1 La #5 5.6 
5 2 1.00 5:3 
6 2 1:25 2-3 
7 2 1.50 4.9 
8 2 1:75 4.8 
Step 4: Run program: 
USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT NO. DATA PTS. 
8 INPUT NO. BASIS FUNCTIONS 
4 INPUT NO. INDEPENDENT VARS. 
2 INPUT F1i(X(1,1),...,X(N,1)) 
1 INPUT F2¢X(1,1),...,X(N,1)) 
EXP -X(1,1) INPUT F3(XC1,12,...,XCN, IDD 
EXP -X(2,1) INPUT FACX(1,12,...,X(N, ТУУ 
EXP (-Х4:1,1)жХ(2,1)) INPUT X1(1) 
а... еменнен Ф-Т рЫ ир с 
1.00 INPUT Y(1) 
6.2 INPUT Х1(2) | 

enter data from table іп step 3 
1.75 INPUT У(8) 
4.8 {wait about 33s} (results printed) 


Step 5 (PROBLEM 2): Fit the data in the table below to a third- 
degree polynomial 


F(X) = Wy * WX * w4X* * w, X? 
I XCIO: XXL. EM 
01 1.02 
05 Ladd 
10 1.23 
2 169 


O <Y O UO FF U Np 
Q N = мі 
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Step 6: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN INPUT NO. DATA PTS. 

8 INPUT NO. BASIS FUNCTIONS 

7% INPUT NO. INDEPENDENT VARS. 
INPUT F1i(X(1,1),..., X(N, I) 
INPUT F2(X(1,1),...,X(N,I)) 


X(1,1) INPUT ESKI T рол CN) 
X(1. D *x2 INPUT FA(X(1,10,..., XCN, Г 
Х(1, Г) жжз INPUT ጂጊ ‹ ጊ ን 
. ()1 INPUT Үс12 


1.02 INPUT ጂጊ ‹2ን 


enter data from table in step 5 


5.2. INPUT Y(&) 
655 {wait about 50s} (results printed) 
دت‎ ег printed) |___ 


Step 7: Discussion: 


Error - 
Та) 
is minimized by REGRP. The choice of basis functions is, of course, 
critical to the ا‎ problem. For example, a regression to 
one basis function, 1, could be done in problems 1 and 2. A 


solution will be нн. E the total error will be large even 
though it has been minimized. Therefore, a choice of basis functions 
should reflect the nature of the data, or an expected result. Since 
the data in problems 1 and 2 was fabricated, the choice of basis 
functions was simply those used to construct the data points. 

The number of data points in both problems was rather small. 
For some regressions, this may be a problem, especially when the 
number of independent variables is large or the number of basis 
functions is large relative to the number of points. 


Step 8: Suggestions: 


Record some data. For example, take your pulse before, during, 
and after some stressful activity, (like a viewing of the 6:00 news). 
It might look like figure 6-6. Choose a number of basis functions 
which reflect the nature of the data, a displaced "hump", centered 


at time t. 


nok dg 
ғы, ENG 
ct oct 
rot 
ct ct 


PULSE RATE 


Use REGRP to find a regression curve, and plot it using PLOT2. 
If these basis functions don't work, choose some other set. You 
may want to regress part of the curve at a time. For example, the 
initial rise is much different from the final fall in figure 6-6. A 
better choice for the basis functions might be а set of exponentials. 


NNNM 
Е 1፳8፳8588588888፪8፪፪.. 
et LE tt Фа 


STRESS ON p^ OFF 


Figure 6-6 Pulse rate response to stress 
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Demonstration 6.5 LINEAR REGRESSION WITH PLOT 


‚2 LET ЖЕЛГЕ 
583 LET MXMUT=Siae 
584 LET РОЇМТа 2257 

р 


C] (ቨ CT] 

AD D D 6) 5፡5) 9) ህ)ርዐር 
BIP БЕЛІГІ Cay 

| PRU 

IAA 


3 33 AICA H T] 
44 


SOOO Cà To = IO CÛ ~, 
ж 
Ini 


chen ca on an сл oua die 
79:8 nro 


ከ1 
ew 
THT 


г- 
г 


ok CHER 


ሂክ ሂክ 
Tag euo а 


тата ж ~ ра 


Step 2: Enter subroutines: VINAX POINT.1 SCALE.2 PLOTD MXINV 
MXMUT REGRP 


Step 3 (PROBLEM): Find the slope and y-intercept of the least- 
squares fitted line for the data in the table below: 


I: KCI): YCIÓ: 
1 2.0 0 

2 3.0 0 

3 3.2 0 

4 4.0 1: 
5 Aad 4.0 
6 4.5 4.5 
7 6.2 6.0 
в 6.5 6.1 
9 6.8 9.1 
10 7.0 10 


Step 4: вип ргодгат: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT NO, PTS. 

10 X(1),Y(1)=? 

2.0 

Q X(2),Y(2)=? 


enter data from table 


9.1 Х(10),Ү(10) =? 
——t 
10 {wait about 15s} (results printed, plotted) 


Step 5: Discussion and suggestions: 


The linear-regression problem, aS opposed to the problems of 
demonstration 6.4, have the two basis functions 


Е = 1 
1 
F5 = K 


The regression curve is a straight line: 


F(X) = Wi * м2хХ 
апа w, із the y-intercept and у, is the Slope. 

Exponential and power seriés regressions may be done simply by 
entering the appropriate basis functions into lines 593 and 594. 
Use different basis functions to regress the data in step 3, 
Compute the error for each set of basis functions, (see spec-sheet 
for REGRP), and see which choice best fits the data. 

Find some real-world data for a linear regression, or an 
exponential regression. The regression line or curve may be used 
to make a prediction, given the past set of data points. 
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EZOO REM REGEN (NLLS REGRESSION? 
S202 DIM Oth 

5204 FOR Ја] ТО M 

Бабе. РЕТНТ "ІНРИТ INITIAL БЕГІС: 
EE 


6208 INPUT uii 


2214 LET üzü-1 


6216 CLS 
2218 PRINT "ITER. "; 0 
8220 PRINT "J: PARI: 


ЕЕгЕрР FOR 151 To አክ 

B224 PRINT Jj TAB 40110) 

БЕРЕ LET Cid (ሀ + ВВЕ HI „т 
Alli ¥. Ui 


"CONTINUE?" 


xm 
ен" . THEN RETURN 


сто TI Yt 
Па pa to hy ut 


Е: 


2252 LET ПІН БЕЙ 
S254 LET GIJ =[ህዉ F%-F ив 


LET Wi =M 


፲ኽ cu cn en cn en ሂካ cu en en ርሽ en pn ei en (en en e Gre ርክ 
ኪን o po ከን ከን йд ከነ PO ha ከነ pcr ኪነ po ከ፲ ከ) ከን ኮን ኪን ከ፲ ha ከጋ 
iD uuu uo e ca ко Ca en -4--4 4-4 4 Chi ТЕТІ Ch ሂበ СЛ 
тре fa C Ch T=- POL C CR. {= ከ) 63 Cü rt p= TO P Ch Ti 


1 11 SUE 

| LET МЕМ 

4 EET Беј 
LET s$s" iI Cl ጅጅ КМ)" 
со SUE MAMUT 
FOR 151 М 
КЕТ чыңа БО вата, 11 
HERT ʻi 
50 Ta 2214 


Figure 6-7 Non-linear least-squares regression subroutine 
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REGRN 


This subroutine performs a non-linear least-squares regression 
of a set of data points 


СОК. Goss ok 


11 Е, 1 ір Ер Р 
То a function containing М parameters о той 
Бог example, F(X,,X,) =w, X, + W. X. +w. X contains 3 
parameters. E ue bot 2 2 3 1 2 і 
The weights w сму are found, such that the error 
Error = 


is minimized. 

An initial guess for the Parameter values is entered, and the 
subroutine iterates to a solution. At each iteration, the current 
parameter values are printed. Convergence is indicated by non- 
changing parameter values. 


CALL: P = Number of data points 
E = Number of independent variables: X... ХЕ 
N = Number of parameters: w per WN 
FS = Regression function: віх ЗЕР ИЙ SS АМ 
Х(,) = Array of independent variable Valües of data points 
Y€) = Array of dependent variable values of data points 
WO = Array to be returned from REGRP with parameter 
values. Must be dimensioned . 
DIM W(N2 
REGRN: INITIALIZATION: 
Prompt: 
"INPUT INITIAL PAR(I)" 


(where I = 1,...,N) 


RETURN: WC) = Regression coefficients 
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Algorithm (non-linear least-squares, thanks to Anant К. Agarwal): 


Given: 
1. P data points 
| ((Х1,Ү,),..., СК Ур?) ; Xi = gE 
| 2. Ап approximation function 
| А 
| э =: 
| FCK;W) ; W = (Ws cco 


3. A current с. В, iteration estimate for the parameter vector: We 


The new estimate, W , is 
L+1 


" m (де ,95፳›... «Фк де 4-1 ү ocn Де, 
Weer "| ዐግ Qi 091 0% 0% 
«де „де Те «Og. ይ » : су Ов > 
дяк 0% ዕዝ OYN Oy 


aera DE „де у = Ò FOD Ж де) 
9% ወ 09%) дяк 


Ora. 


gu 


<Y-F „де > = (Y; - Ес.) ] ж 


I?! 


————— መታ ፍው 
Q 
= 
с 


| 
n 
г 
^ 
% 
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Demonstration 6.6 GENERAL NON-LINEAR LEAST-SQUARES REGRESSION 


Enter main program: 


ሠ 
ፎኑ 
(D 
ке) 
፦ 


818 REM рЕМОВ.5 
211 LET REGRH=5200 
512 LET MX INL 55088 
Es БЕ МАМИТ авіда 
514 INT “INPUT MO. DATA Ars" 
545 INPUT р i 
. РІВ PRINT "INPiiT НО. PARAMETERS 
617 І 
518 P 
T VARI 
519 Th 
626 г 
521 D 
622 p 
ваз P 
б „11; 
624 T 
625 ғ 
626 Ci 
БЕТ E; 
| PRINT Li ig 
3 INPUT I | 
та 


Step 2: Enter subroutines: MXINV MXMUT REGRN 


Step 3 (PROBLEM 1): Fit the data in the table below to the function 


F(X) = (1 + w Х) 2 
І: X(I): YI 
1 . 20 1.02 
2 .40 1.05 
3 .60 1.07 
4 . 80 1.09 
5 1.0 1.11 
6 1.. 2 1.12 
7 1.4 1.14 
8 1.6 1:26 
assume initially Wy = “> = 1 
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Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT КО, DATA PTS. 

8 INPUT NO. PARAMETERS 

2 INPUT NO. INDEPENDENT VARIABLES 

1 1МРӨТ_Е‹Х‹(1‚,1),...,Х(Е,1);%#(1›,... (муу 
(1+W(1)*X(1,1))**W(2) INPUT X1(1) dx neam 
72 INPUT Y(1) 

1.02 INPUT X1(2) ———R—M HÀ 


enter data from table in step 3 


1.6 INPUT Y(8) 
1.16 INPUT INITIAL РАВ . . — ፡ — 7 
1 INPUT INITIAL PAR(2) 
1 (iteration O printed) 
CONTINUE? 
(enter) (iteration 1 printed) 
CONTINUE? 
keep iterating! 
{enter} (iteration 8 printed) 
CONTINUE? 


і N (program stops) 
| Step 5 (PROBLEM 2): Fit the data in the table below бо the function 
FO,,X4 = w, (X, - Wo) Ха 


1 + We Xi * W, Хо 


፦ 
9 


.ЗЕ-6 
.6Е-6 
.9E-6 
.ЗЕ-5 
DES 
.OE-6 
. 6E-5 


HON YO OF WON) н 
Oro Oc OO NN N N N| K 
O O O O O Ul UO N Ol O 
ав вен 

N OY UI N N O x ON] ж 
QI O ‹ቨ O Q Oi O OT O ጩ 

Q ON Hp = = OQ O| < 


0 


Assume initially: w,-.01,w 5.1,445.1,М,:.01 
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Step 6: Run program: 


USER ENTERS: COMPUTER RESPONDS: ; 
RUN INPUT NO. DATA PTS. «Й 
10 INPUT МО. PARAMETERS | 
& INPUT МО. INDEPENDENT VARIABLES — | 
2 INPUT FCK(1,1),...,X(E,1);W(1),... WIND 


rn መ መ ЛА 
W(1)*(X(1,1)-W(2))*X(2,I)/(1+W(3)*X(1,I1I)+W(4&)xX(2,I)) 
INPUT X1(1> 1 


2.5 INPUT Х2 (1) 
ረጋ ን INPUT Y(1) 
3.3E-6 INPUT Х1(2) 


enter data from table in step 5 


12.5 INPUT Y(10) 
БЕС ең) айы INPUT INITIAL PAR(1) 
. 01 | INPUT INITIAL PAR(2) 
¿l INPUT INITIAL PAR(3) 
s.l INPUT INITIAL РАВ(4.) በጠ 
‚01 (iteration б printed) 1 
? 
{ ENTER } (iteration 1 printed) 
CONTINUE? 


keep iterating! 


tenter) (iteration 7 printed) 
CONTINUE? 


N (program stops) 


Step 7: Discussion and suggestions: 


The non-linear least squares problem usually requires many 
iterations before a solution is found. There are many cases when 
convergence does not оссаг or when the parameters get set to 
unrealistic values. That is, it is trickier to use REGRN than 
REGRP. On the other hand, if a theoretical formula is to be 
applied to a set of data points and approximate values of the 
fitting parameters are already known, then REGRN may give good 
"optimized" values for the parameters. In this problem, as in 
the root-finding problem solved by ROOTN, the choice of initial 
values for the parameters is very critical: the results from one 
initial guess may differ radically from the results obtained 
Starting with another initial guess. 

The demonstration might be improved by limiting the values 
of the calculated parameters. For example, if the "realistic" 
values that might be taken on by a parameter м, are between 0 and 2 
then the subroutine would limit the value of “| to between 0 апа 2 
at each iteration. 

Use some different data, plot it using PLOT2, (Chapter 2), and 
try to fit it with various functions using REGRN. 
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Demonstration 6.7 EXPONENTIAL. REGRESSION WITH PLOT 


Step 1 (PROBLEM): Given the data in the table below, find the 
parameters Wd and Wo of the least-squares fit to the function 


F(X) = м EXP (w, X) 


1 
Таз X(T): УС 
1 SQ 2 
і 2 22 4 
1 3 135: 6 
4 «55 12 
5 .76 20 
6 „81 23 
7 .96 36 
8 .98 38 
Step 2 Enter main program: 
вай БЕН Di 
541 LET Я 
242 LET H 
543 LET H 
4 LET P: 
Е LET Fi 
вав LET wu 
547 БЕТ 
242 n 
I 549 LET 
кәр LET 
551 БІН 
252 ЏЕТ хр 


ПЕЙ 
148 ቫቱ Ko 


CT Да በ1 በበ ጠ 


IO ri ~ 


ብ FP 
mn. алыш 


Hi Ki — +: — (Tt 


fap: E зо En 


Step 3: Enter subroutines: VINAX SCALE PLOTD POINT MXINV MXMUT 
REGRN 
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Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT NO. PTS. 

в X(1),YC1)=? 

.01 


2 Х(22,Ү(25-? 
A д, 2, 


enter data from table іп step 1 


36 X(8),Y(8)=? 


.98 

38 INPUT INITIAL PAR(1) 

ыз INPUT INITIAL PAR(2) 

> K ыр- ر‎ 
(iteration 0 printed) 


2.5 
هه‎ CONTINUE? 
{enter} (iteration 1 printed) 
CONTINUE? 
keep iterating! 
Center} (iteration 5 printed) 


CONTINUE? 


CONIINUE? | (20‏ سس 
N (program згорз)‏ 


Step 5: Discussion and suggestions: 


This demonstration makes use of 7 subroutines, making for a 
short main program. Since the form of the fitting function is 
predetermined, fewer prompts are required than demonstration 6.6. 
After the parameters are found, а plot of the fitting function is 
done, along with the input points. 

Try different initial guesses to determine i£ more than one 
set of "optimal" parameters exists. 
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Chapter 7 DIFFERENTIAL EQUATIONS 


We all possess some ability to predict the future. When we 
do this, we synthesize our Knowledge of the state of the universe 
at the present time with our knowledge of how it is likely to 
change, or how it has already been changing in the recent past. 
Thus, we solve a sort of boundary value problem when we perform 
such a prediction. 

Differential equations are encountered in problems dealing 
with the change of a variable or variables in time or Space, or 
with respect to some other set ОҒ parameters. Тһе solution to 
such a problem is called a "trajectory": the predicted "path" in 
space, time, or some other parameter, taken by the variable. The 
boundary value problem consists of: 


1. Knowledge of the present state of the "universe": The initial 
conditions. That is, the value of the variable at t=O or 
Х-О, or 


2. Knowledge of the way the "universe" has been changing: The 
differential-equation(s) describing the mechanics of the 
problem. 


The subroutine RKUTT solves a system of ge order differential 
equations, with initial conditions, using the 4 order Runge-Kutta 
algorithm.  RKUTT can deal with a vector of output variables, 
changing with respect to one input parameter. Differential 
equations involving, higher ordered derivatives may be resolved 
into a system of 1 order differential equations. An example is 
given in demonstration 7.1, problem 2. 


TI ФК СД, К-т та“ 


“መ፲ ትኪ 4:61 даа 


` E) 6.63 C3 б) 53 ርሀ бу = 63ር 6 2 63 6) 
келіс ሆላ ጄነ ርነ ከነ ከነ ከነ ኪንከንዮ። s ፦፦ 
B) ዉቶክ;6 00 си P NJ ба Gj Tı {= [О G 


ول“ ل ча dad‏ 04 ل ل 74747474 


А 
7248 
l BE 
40 МЕД, Беда altel, ІЗ ЖІКІМІ,27 

+E (dl. 41) ИВ [FE J Bik [ፌ፤, 21) 1 #83 
TELE NEXT J 
7044 НЕХТ І 
T8456 RETLRH 
7048 FOR Jal Ta E 
FES LET Wt svAL =ዌ 
FASE НЕХТ 4 
72854 FOR Ј=1 ТО E 
rese LET Kig,K)zHsUmL Fe tu 
това МЕХТ J 
TOGO LET К=К+1 
?ӘБЕ RETURN 


Figure 7-1 Differential equation solving subroutine 
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RKUTT 


This subroutine numerically solves first-order vector differ- 
ential equations of the form 


ау, OO РОМ OD, ... YN CKD) 
ах 
ау CX? Е OG У, OO... У X0) 
ах 


over the interval: Х,<-Х<-Х The values of the solution vector 
are specified at the Tower ейа- -point of the interval: 


Y, X) > а, 


Ук FR? ay 


The subroutine uses a fourth-order Runge-Kutta algorithm. 
Values for the solution vector are calculated at each step in the 


independent variable 


= CI-1)*H ; = Ses - . 
Xi Xp + )х 1 125 A INTC (Xy X; 2/8) + 1 


where H is the user-supplied step-size for the integration. 


CALL: E = Number of variables: (Y) 
- Number of equations: {Үү'= F.) 
= Lower end-point: X 
= Upper end-point: X 
Е ሀ 
= Step-size 
ЕФ() = Array of functions: (FU 
YO) = Values of ҮІ ат Хех, 


RKUTT: Changes I J K K(,) М У(,» X X$ 


T шр 


RETURN: N Number of steps 
ሃ(,› = Solution at all steps: 


МСЈ,1» = ¥ ሻኻ + СІ-1)хжН) 
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ዴር ሽር 


Algorithm ca th огаег Runge-Kutta): 
Given the solution point.at a particular step 
(X;Y ,... YE) 


the calculated solution point at the next step is 


ረ ops. АҚЫ ж e 


where 
Y; = 1/6 (кұс) + 2крсј) + 2፪,(7) + к.с) 
Ki) = Ах ВО бр 
KJ) = Ах F OC AXSY HE (1)... УК CE) 
K4) = Ax FOGG AX; Y BE CO, LLL У А (ED) 
K,GD = Ax вус Axiv,- Ка ም ከን KED 
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Demonstration 7.1 PLOT COMPUTED SOLUTIONS TO SEVERAL 


SYSTEMS ОЕ DIFFERENTIAL EQUATIONS 


Step 1: Enter main program: 


670 RE 

671 LE 

672 PR EQUATIONS" 
673 IN | 

274 FR LOWER, UPPER L 
ІМІТЕ" 

675 

576 

577 

578 

679 

вав 

851 

582 % 


Step 2: Enter subroutines:  RKUTT 


Step 3 (PROBLEM 1): Plot the solution to the second-order 


differential equations: 
3 
dY (X) -5Үү+ү У 0) 1 
ах = 
dY,(X) C EXP(2X)xY, Y, (02 1 
ах 
over the interval: О<=Х<=1 


ዘ 


System of 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN INPUT NO. EQNS. 

2 INPUT LOWER, UPPER LIMITS 
О 

1. " INPUT STEP-SIZE 

202 INPUT БІ(Х,У(1),...,У«(М)» 
-ዱቾኛ‹ርጊን+ኛ‹(መን›።።ጋ INPUT Y(1) АТ Х-0 

1 INPUT F2(X,Y(1),...,Y(N)) 
-EXP (2жХ)жу (1) INPUT Y(2) AT Х-0 

1 (wait about 1ቪ503) (results plotted) 


Step 5 (PROBLEM 2): Plot the computed solution to Bessel's equation 


X° d*y(X) + X ду (СХ) + (X? - му = O 
ах? ах 


for М = 0, over the interval ዐ‹ጂ‹=5 
То ао this, let 


ҮҮ = Y 
Y, = ах 

2 ах 

This gives the system of 155 огаег differential equations 

Die S 

аа ах 2 

аү 

2 mo e > ах 2 ма 22 = = 

ах = dX? = ፦ | ах + (X Мао /፳ = -Yo/X Yi 


У, is the solution to the original differential equation. Using 
the initial,gonditions Y (0) = 1 ; Y4,(0) = O, the solution will 
be the zero -order 865561 '5 function of the first kind, Ј_, (see 
chapter 3). Also, since J.' = > J,, Y, will be the inverted js -order 
Bessel's function of the first kind. 2 


Step 6: Run program: 


USER_ENTERS: COMPUTER RESPONDS: 
RUN INPUT NO. EQNS. = 
2 INPUT_LOWER, UPPER LIMITS 
ааа INPUT LOWER, UPPER LIMITS ` 
1Е-7 


5 INPUT STEP-SIZE астананың 
05 INPUT РІСХ,УС1),...,УСМ)) 


ee РАКА, СІ), ሬር YONI) | 
INPUT Y(1) АТ Х-1Е-7 
کے‎ EU УФІ) АТ ХЕ _____ 
1 : INPUT F2(X,Y(1),...,YCND) 
— <<< ر‎ ee m 
-Y(2)/X-Y(O1) INPUT Y(2) AT X-1E-7 


——F— .— h|A[ YU Á‏ > ] ددد 
О {wait about 2203) (results plotted)‏ 
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Step 7: Discussion and suggestions: 


Problem 1 and problem 2 were both second-order problems. 
Problem 2 showed how a differential equation involving both first 
and second order derivatives may be resolved into two equations 
involving only the first derivative. Differential equations 
involving higher order derivatives may also be resolved into a 
system of first-order differential equations in this way. 

There are a number of algorithms which are used to numerically 
solve differential equations. The Runge-Kutta method integrates the 


functions F_,...,F, at each step in the independent variable, X, 
using a rulé similar to Simpson's rule, (See chapter 4), starting 
from the initial values of the dependent variables Y,.,...,Y.. 

For example, the 1 order problem is: 1 М 

ау сю _ , қ x: 

ах . - = F(X,Y) . YCA) = S 
Integrating: | А+АХ Аз? АХ x 

ЧОО = %,» F(X, бай + JFR, аў +... + JFR, Yak 

А А+ Ах Avtar 


where ¥ is the value of Y in the small interval. Each integral is 
calculated using a rule. The integrations are continued in this way 
until X reaches its upper limit, B. - 


Use RKUTT to solve some other differential equations. Try some 
higher order systems of differential equations, and some differ- 
ential equations involving higher orders of derivatives. 

«Тһе general resolution of ап М" -order problem into а system 
of 1" -order differential equations is: 


1. original equation: 


I 
) (,CO 8-00 „кх „о 
ах 
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2. Тһе М first-order equations аге: 


uns 128 

s жу. 

dX 

Š N 

ау, _ _ | FOO у du YN 
ах Gy OO ден Gu CO 


* Р th = у 
The initial conditions are specified, by the N^ -order derivatives 


of the function Y at A: ҮКА? = ү (А? 


Demonstration 7.2 PLOT COMPUTED SOLUTIONS ТО SEVERAL 


FIRST-ORDER DIFFERENTIAL EQUATIONS 


Step 1: Enter main program: 


Step 2: Enter subroutines: RKUTT VINAX SCALE.1 PLOTD 
Step 3 (PROBLEM 1): Plot the computed solution to 


ЧҮ = -XY $  Y(-5) = EXP(-i2,55 + -5<-Х<-5 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 


UN INPUT INTEGRAND-F(X,Y(1)0) 
-ХжҮ<1) INPUT LOWER, UPPER LIMITS 
35 


I Х--5 


ЕХР -12.5 {wait about 505) (results plotted) 
———  ШЗЕКЖ መመር 305) (results plotted) 0 


Step 5 (PROBLEM 2): Plot the computed solution to 


dY = -XY* 


= » 


а 


УСІ) & 2 ; ጊ‹=ጂ‹=2 


ж 
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Step 6: рип ргоагат: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT INTEGRAND=F(X,Y(1)) 
-X*YC1)**2 INPUT LOWER, UPPER LIMITS 
1 
2, INPUT У(1) AT Х=1 
2 {wait ароці 1т155) (results plotted) 
Step 7 (PROBLEM 3): Plot the computed solution to 
ам = -Х/Ү н ፕ(0) = -1 2 ዐ‹=ጂ‹=1 
ах 


Step 8: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN INPUT INTEGRAND=F(X, У (1)) 
-K/Y (1) INPUT LOWER, UPPER LIMITS 
о 

1 1 1 = 

-1 (wait about 50s} (results plotted) 


Step 9: Discussion and suggestions: 


The last demonstration main program could also have been used 
to plot out the computed solutions to these problems. However, 
less prompts are required for this demonstration, since the order 
of the system of equations is pre-determined. 

The problems in this demonstration may be solved by simple 
analvtical methods, such as separation of variables. Compare the 
numerical results with the analytical results. Use different step- 
sizes, and observe the effect of varying the step-size оп the 


accuracy of the computed results. 


Спаріег 8 NUMBER REPRESENTATIONS 


There are many ways to symbolically represent a number. 
The ZX81 and TS2068 use this internal representation for "floating- 
point" numbers: 


8 bits 1 bit 31,bits 
А, 
ЕХРОМЕМТ SIGN FRACTION 
where а "bit" is either 0 or а 1. То get back to scientific 


notation, assuming EXPONENT, SIGN, and FRACTION are all integers 
in base 2 notation: 


number = (-1) (SIGN, 2 (FRACTION, „#2 72+27 ^) ко! (EXPONENT, о-128) 


where the subscripts "10" indicate base 10 notation corresponding 
to the base 2 integers, 

This chapter presents a few other ways to represent numbers, 
including ratios of integers, 9's complement notation, and different 
base notations. 

RATIO iS a program that finds two integers whose ratio is 
close to or equal to the value of а given number. The real number 
system is composed of both rational andirrational numbers. The 
rational numbers are those whose values are given exactly by a 
ratio of integers. Since all numbers are represented by the 
2Х81 and TS2068 with only up to nine digits (base 10), all of 
the numbers which are represented on the computer are essentially 
rational: for example, 


.123456789 = 123,456,789/1,000,000,000 


ВАТТО, помемег, looks for the smallest integers whose ratio із 
within 10 of the given number. 

The subroutine OPERA extends the accuracy of the computer for 
multiplications, additions, and subtractions. Operations with 
very large decimal-place accuracy take much longer than the usual 
9-digit accuracy operations, (since the program is in BASIC, rather 
than in assembly language). Тһе subroutine uses a similar rep- 
resentation of each number to the ZX81/TS2068 internal representation. 
However, base 10 notation is used rather than base 2. OPERA is 
used in demonstration 8.2, a reverse- polish notation calculator, 
with user-defined precision. 

The program BASER changes the representation of a number in one 
base to another base. Both integers and fractions are handled by 
BASER, and fractions are truncated to the same accuracy as the 
original representation. . 

Finally, ZXCAL turns the ZX81/TS2068 into a scientific notation 
calculator, with successive results scrolled onto the screen. 
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2285. РЕН PROGRAM RATIO (RATIO-FI 
2255 PRINT “INPUT HUMBER "; 


Qo co en 
ооо 
ноз 
Gon 
гост 
и n 
H 
= 
wi- 
ин 
Tir то 


8012 LET сав 
25014 IF ABS (2-Е! 16-86 THEM GO T 


24 IF С<Е THEN LET ፳=ጸ+3 
26 IF С›Е THEM LET В=Ё+1 
28 GO TO за15 


Figure 8-1  Ratio-finding program 
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RATIO 


This program finds a ratio of two positive integers which 


approximates а user-supplied positive real number: 


- E <= AD. 9 


A 
B 


where A and B are integers, and E is the user-supplied real 


RUN: 
RATIO: Prompt: 

"INPUT NUMBER" 
STOP: - Input number 


Numerator 
Denominator 


E 
А 
B 


ዘ 


E ል B and А/В printed 


number. 
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Demonstration 8.1 FIND RATIOS FOR SEVERAL IRRATIONAL NUMBERS 
ACCURATE TO 6 DECIMAL PLACES 


Step 1: Enter program: RATIO 
Step 2 (PROBLEM 1): Find a ratio for 4. 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
PI {wait about 45s} (results printed) 


Step 4 (PROBLEM 2): Find a ratio for e. 


Step 5: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
GOTO 8000 NPUT IMBER 
EXP 1 {wait about 5m45s} (results printed) 


Step 6 (PROBLEM 3): Find a ratio for 43/2. 


Step 7: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
8000 N 
SQR 3/2 {wait about 3m15s} (results printed) 


Step 8: Discussion and suggestions: 


It's interesting that the ratio of two fairly small integers 
gives a close approximation to the value of . It takes about 7 
times longer to find a ratio for e. There are, of course, many 
more irrational numbers which may be subjected to RATIO. In fact, 
the set of real numbers might be described as the set of irrational 
numbers with a few rational ones thrown in. However, as these 
examples show, given any irrational number, a rational number may 
be found which is arbitrarily close to it. This is fortunate, since 
computers cannot yet grasp the concept of the irrational number. 

Try RATIO on an integer, a number with a repeating digit, such as 


.0033333333... ог .457457457... 


and other numbers which may ог may not be rational. 
Find ratios with less accuracy than 10ppm by changing line 
8014 of the program RATIO. Print out all results which are within 


a certain range of accuracy. 


broutine 


nm ЏИ] 
lu ፦ m 
й. [1 са 
[шї] A па 
in ul 

" ይና = б 
га HI! ዘ É 
Iu ЈЕ T ሂጋ 
፲፪ к ре 
fale - І 

> lii [መጁን 
- і 24) md 
к 2 1F 
m ዘ ዴታ 
Бај Р ሀ] 
= TE 
БЕ HE 
т ር 
її + 
ІН 
п. ~ Г FZ W H onn 
let ° нт "ETE | uc ud ee Ш ው የ] 

3 т 1 к H H ኤ . | 
፳ ሁሁ ኑዮ፦ዬጅዬ db BO pO መዚወሁኑ'ሩ LE XEL врх ወ. итар КЕЕ ЕЕ 
ii ШЕШШ Шш Dx 


Ш — ЛА НО BA bin. wc Bret ICL. РР РУ Р А РР РУ HC същ ши ПР JHT E ubi 
للك‎ JAJA SA clk нина GOLA ZE AAA C mI EE ENE ПЕШО р 


к 
Ten ዘ1 O HD Fade шенге TA CJ 1) በ፲) ርሀ ርቄ C <f iD O GÛ АТ е АУ ВЕСИ а а DE са С <f 
OHA C EA A A A КАТ 0| reri cerne ቁ CIF a FF <F IZ Ii Û lı tN IG IO i ID с P= f EU A 3 е A A A e Aa CU Cd O 
Яна са c C ed e 4 eT ና45፡1ና1ና1፡15151፦ да E г р 1 
YO E а ату] 03 03 mr ከ1 ccv mc ro a a OD IC 030 


1 0 NO O x 3) 03. 0 


Figure 8-2 Multiple-precision operations su 
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F+ p> ра ра ре бу C) Ca Ч Са З 
f £f F [U бы Cü Cf P. — D С) 


"3 C 


1 £r са Ca 01 Di Ға 01.0) ca po ca ርዕ ርዕ ር) Сі) 


CH GI са Ca CI CO Ca C3 3 


ር) 


га са са £0 “J TE f 4 


[x 
es 
53 
A+ 
2378 
2350 
8388 
‘8384 
eng 
2402 
2484. 
8406 
2488 
2418 
8412 
8414 
8416 
8418 
a4zü 
54== 
2528 
2502 
2504. 
a506 
esas 
5535 
5512 
2514 


"nr 


H 


E 
ATT +H 


Tuo 
moe mto 
ло 444 


ha 
ж] 


{ 
4 


= 
um, 


ኀ 
“ገ 
ШЦ TY 


т 
” 
~ 


«Шан 


кашы ее тү түш 


1. 
ж 
ы 


-ፒ፻ጋዘ፦ በ፤ ዘ4ሥ ዘ 


H 
2 
1 
1 
T 
I 
Ih 
I, 


PRINT "4%; 

FoR Ісі ТО в 
PRIMT Z(A-I+1,1:; 
МЕХТ І 

PRINT "E";zin-ez,.il 
RETURN 


ОРЕВА 


This subroutine performs multiple-precision multiplication, 
addition and subtraction of floating-point numbers. The precision 
is caller-defined. The subroutine is designed to be used іп а 
Reverse-Polish Notation (RPN) calculator, as shown in demonstration 
8.2. 


CALL: М = Precision: number of decimal places = Nx8 
Z¢,) = Register array. Must be dimensioned 
DIM Z(8*N+2, 2) 


X$ = Number in scientific notation 
УФ = operation: 
"у" for ENTER : Put number in accumulator 
"+" for addition : Add number to accumulator 
"-" for subtraction : Subtract number from 
accumulator 
"ж" for multiplication : Multiplv number by 
accumulator 
OPERA: Changes А AS B C D E F G H I J XO YO 26,2 
RETURN: Operation performed. Result is put in accumulator. 
Accumulator printed. 
Summary of subroutine operation: 
The register array has this format: 
Register name: Variable: Format: 
Accumulator 2‹።.,1) [Е ЕЕ... ЕЕ Е ЗЕ | 
х Z(x,2) LEEF... F F Е[з]Е | 
8N places 


Where .FFF...FFF is the fractional part of the number, S is the sign 
of the number (+1) and E is the exponent of the number. The frac- 
tional part of the number is adjusted so the first number after the 
decimal point is non-zero. 


Multiplication (MULTX): 


Lines 8200 to 8214: Multiply all digits of accumulator with all 
digits of X and store result in YO. 

Lines 8216 to 8226: Transfer portion of У() that is greater than 9 
to next digit, etc. 

Lines 8228 to 8234: Look for first non-zero digit in YQ). 

Lines 8236 to 8244: Transfer Y() to accumulator such that first 
digit after decimal-point is non-zero. 

Lines 8236 to 8248: Adjust Sign and exponent. 
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Addition/Subtraction (ADSUB): 


Line 8302: If operation’ is subtraction, then change sign of X 

and proceed with addition. 

Lines 8304 to 8310: Depending upon difference in exponents of X and 
accumulator, determine the register to be shifted (H), апа by how 
many places (Е). Call shifting subroutine, (SHIFT). 

Lines 8312-8324: Convert to ten's complement if sign of accumulator 
or X is negative. 

Lines 8326 to 8336: Add X to accumulator. 

Lines 8338 to 8342: Ten's complement accumulator if result is 
negative. 

Lines 8344 to 8354: Left-shift accumulator to give a non-zero digit 
after the decimal point. 


Line 8356: Adjust sign: 
ном - MN and "а" E мп то "ел! = кы апа "-1” = во” 


Demonstration 8.2 REVERSE-POLISH CALCULATOR 


зтер 1: Enter main program: 


Step 2: Enter subroutines: OPERA 


Step 3 (PROBLEM 1): Find the number of radians in one degree, 
24 decimal places. 


radians/degree = 4/180 = Я ። 0.00555 


Step 4: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
RUN INPUT PRECISION 
З ІМРОТ МОМВЕВ 
3.141592653589793238462643 INPUT OPERATION 
> (number is printed) 
INPUT NUMBER 
.0055555555555555555555555 INPUT OPERATION 
* (result is printed) 
Step 5 (PROBLEM 2): Using 15 place accuracy, compute the sum 
S= 1+ X+ X? кк + х* + xX = X(X(X(X(R+1)+1)+1)+1)+1 


for X - .005 


to 
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Step 6: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

RUN INPUT PRECISION 

2 INPUT NUMBER 

1 INPUT OPERATION 

> (accumulator is printed) 


INPUT NUMBER. ሎሎ 
INPUT OPERATION 


005 РИ OPERATION, aa 
* (accumulator is printed) 


INPUT NUMBER 
1 — == | INPUT ORRERA LILON aa 
Т (accumulator is printed) 


INPUT NUMBER 


————— H h 


4005 ТКРИТ ОРЕВАТІОМ 
ж (accumulator is printed) 


INPUT NUMBER 


E E‏ ا 
INPUT_OPERATION‏ 


፦ 


+ (accumulator is printed) 
INPUT NUMBER 

а NUT OPERATION, ee 

ж (accumulator is printed) 


INPUT NUMBER 
1 INPUT OPERATION 
+ (accumulator is printed) 
INPUT NUMBER 
.005 INPUT OPERATION 
* (accumulator is printed) 


INPUT NUMBER 


1 ІМРОТ ОРЕВАТІОМ 
+ (accumulator is printed) 


INPUT NUMBER 


.005 INPUT OPERATION 

ж (accumulator is printed) 
INPUT NUMBER 

1 INPUT OPERATION 

ра (accumulator is printed) 


Step 7: Discussion and suggestions 


Problem 1 and problem 2 keep all of the interesting digits 
that are usually lost in computer calculations. Many multiple 
precision calculations are amazingly symmetrical; for example, those 
involving 9's. 

Appendices B and C give some physical and mathematical constants 
and conversion factors with many decimal places. (The value for 
« has been calculated to many more pages of decimal places than is 
shown in appendix B.) Many calculations can result in the amplifi- 
cation of round-off errors, and in these calculations, the larger 
the number of decimal places, the better. 

Try writing a division subroutine for OPERA. Of course, multiple 
subtractions might be done, but this would take much time to perform. 
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; 
BASE "SB 

T 

+ 

x 

Kis". THEN 85 T 


p 

Dn ве 

2524 МЕХТ E 

3625 LET NzK-1 

8628 IF Мао THEN Go То 
8630 CIM Bin: 

S632 FOR ፲=3 To K-41 
8534 (Е T 
45: Е 
Бік = 
KLE N 
з 


литий б 


ew ፓኑ 
Bes n 
Sea 5] 
3555 ሠ 
3558 5 
2578 R 
8672 FOR I 
S674 LET Uu 
S676 LET B 
5575 LET В 
5552 LET 5 
8582 МЕХ 
3584 LET 
5585 IF : 
5555 For 
5528 LET 
RINT CHR 
USE +++ 
8692 NEXT 
8544 RETU 
аура REH 
2702 DIH 


TA и Ген г: 
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Figure 8-3 Base changing program 


BASER 


This program converts a number in one base notation to another 
base notation. The number may contain an integer part and a 
fraction part. Fractions are rounded off to the equivalent number 
of places in the new base notation to the old base notation. 


RUN: 
BASER: Prompts, in order of appearance: 
1. "INPUT BASE(TO)" 
2. "INPUT BASECFROM)" 
3. “INPUT NUMBER IN ВАЗЕ В" 
(where В is the old Базе) 
STOP: Number printed old base and in new base 


Example of algorithm: 


1. Convert 12.3419 to base 2 notation. 
2. Integer part: 


Let 12 = a2? + ва! + с2= + 427 + 


Successive divisions of the integer part by 2 gives: 


6 —D2 о зла de S gies ; remainder a = 0 
0 1. ; 

3 = C2 + d2 я ; remainder b = O 

1 + а29 + 5 remainder с = 1 

О = аж ; гетајпдег а = 1 


Hence, successive divisions by the new base gives remainders equal 
to the digits of the number in the new base, starting with the 
least-significant digit. 

3. Fraction part: 


Let .34 - a2 > b2 ^ + c2 + በ2 “+ 


Successive multiplications of the fraction part by 2 gives: 
1 == == 


0.68 = а + b2 ` + са + 42 жо... 3 integer part а з О 
1.36 5-11 8d + a2"2 +... : integer part b = 1 
0.72 = c + dg ^ +... ; іпседег part c = 0 


ዘ 
m 


; integer part d 


0 1100.0101, 


1.44 = а + 
4. From parts 


N 
0 . 
5. 
о. 
о 
= 
to 
о 
+ 
m 

u 
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Demonstration 8.3 CONVERT SEVERAL NUMBERS FROM 
ONE BASE TO ANOTHER BASE 


Step. 1: Enter program: BASER 


Step 2 (PROBLEM 1): Convert Е.А10С,2 to base 10 notation. 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 8600 INPUT BASE(TO) 

10 INPUT BASE (FROM) 

16 INPUT NUMBER IN BASE 16 

Е.А10С (results printed) 

Step 4 (PROBLEM 2): ‘Convert 10110.00111, to base 10 notation. 


Step 5: Run program: 


USER ENTERS: COMPUTER RESPONDS: 
GOTO 8600 INPUT BASE(TO) 

10 INPUT BASE( FROM) 

2 | INPUT NUMBER IN BASE 2 
10110.00111 (results printed) 


Step 6 (PROBLEM 3): Find the base 16 representation of 1. 


Step 7: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 8600 INPUT BASE(TO) 

16 INPUT BASE( FROM) 

10 INPUT NUMBER IN BASE 10 
3.141592638979323846 (results printed) 


Step 8 (PROBLEM 4): Find the base 2 representation of ў. 


Step 9: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 8600 INPUT BASE(CTO) 

2 INPUT BASE (FROM) 

10 INPUT NUMBER IN BASE 10 
3.141592638979323846 (results printed) 
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Step 10: Discussion апа suggestions: 


Most of us are fluent only with base 10 representations of 
numbers. However, there are at least three other commonly used 
number representations: base 2 (Binary), base 8 (Octal), and 
base 16 (Hexadecimal). These bases are all powers of 2, and since 
all information in a digital computer exists, at least on some level, 
in a binary form, the binary, octal, and hexadecimal number systems 
are convenient representations to use. To change base notation 
within these powers-of-two representations is simple. For example, 
the change from Базе 2 to base 16: 


base 2 number : ... 1011 0110 . 0101 0101 
base 16 number : нз в 6 š 5 5 


the algorithm is to group four of the base 2 digits together, without 
overlapping the radix-point, and writing the base 16 digit corresponding 
to each group. 

Note that an appropriate number of places is chosen in the new 
base notation Бу BASER. This corresponds to the same accuracv as the 
representation of the number іп the old разе. 

It is highly recommended to view the representation of PI іп 
base 32. Note that letters take the place of numbers after vou reach 
base 9, (any base larger than 10). If you get higher than 2, then 
other symbols corresponding to character codes will be used. 
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200 REM PROGRAM ZXCRHL [2981 CAL 
CULATOR: 

82 LET ква" " 

114 REM ses FOR 2:81, USE sss 5 
Шш. 
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= 
EL 
E 
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Figure 8-4 Scientific calculator 
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ZXCAL 


This program makes the ZX81 or TS2068 emulate a scientific 
calculator with the functions: 


+ - ж / LN EXP SIN COS TAN ACS ASN АТМ ЗОР 


The "ENTER" key is equivalent to the "=" key on a Scientific 
calculator.  Successive results are scrolled onto the display. 


Demonstration 8.4 SCIENTIFIC CALCULATOR 


Step i: Enter program: ZXCAL 


Step 2 (PROBLEM): Use ZXCAL to calculate 
EXP( (355/113)*LN(3.2+6.8)) 


Step 3: Run program: 


USER ENTERS: COMPUTER RESPONDS: 

GOTO 8800 (waits for input) | 
3.2 (enter) 3.2 

+6.8 (enter) 10 

LN (enter) 2.3025851 

*355/113 {enter} 7.2337851 

ЕХР (ептег) 1385.4567 


Step 4: Discussion and suggestions: 


It is hoped that this book has shown that the computer is 
orders of magnitude more powerful than a hand-calculator. Hand- 
calculators are appreciated for their portability, but the computer 
can always take over when necessary, аз ZXCAL shows. 

ZXCAL may be included in other programs for a convenient 
calculator for stray calculations. If you have a printer, modify 
ZXCAL so that it prints out the successive results. 
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аррепдїх А ONE-LINERS 


Table A-1 lists some mathematical functions not found on the 
2Х81 апа TS2068, but that may be implemented in one program line, 
as shown in the "implement" column. The center column gives a 
range of values for the arguments, for those who like to stretch 
the limits of their machines. 


Notes: 


1. ТЕ the intent is: Ү-тіп(Х,Ү), then this function may be more 
efficiently implemented by 


IF X<Y THEN LET -Y=X 


2. The MOD function gives the value of the argument in a particular 
MOD. For example, counting in MOD 5 goes Jike this: 


01234012340123 


Table А-1 One-liners 


и = 4х10799 N = 1х1038 
ІМТЕМТ: CAREFUL! IMPLEMENT : 
Z=log 00 BOX«N LET Z=LN X/LN 10 
Z-cosh(X) -88<Х<88 LET Z=.5*(EXP X«EXP -X) 
Z-sinh(X) -88<Х<88 LET Z-.5x(EXP X-EXP -Х) 
Z-tanh(X) -М/10<Х<44 LET Z-(EXP (2*X)-1)/(EXP (2*X)+1) 
Z=cosh | (x) 1‹=፪‹ዝ N LET 2-10 (X+SQR (፪=፪-1)) 
Z-sinh 1сх) -307 Теже N LET Z=LN (X+SQR (K*X+1)) 
Z-tanh +(x) -10X«1 LET Z=.5*LN ((1+X)/(1-X)) 
(0) 2=та псх, v) LET Z=(Y<X)*Y+(Y>=X)*X 
@ г=тахсх, Y) LET Z=(Y>X)*Y+(Y<=X)*X 
@ Z-mod, (X) X,Y integers LET Z-X-INT (X/XYOxY 
Z-remainder(X/Y) LET 2-Х-ІМТ (X/YOxY 
Z-db conversion(x) го H«X«N LET Z=20*LN X/LN 10 
Z-deg to rad(X) LET Z=X*PI/180 
Z-rad to deg(X) LET Z-X*180/PI 
тех“ X may Бе <0 LET Z=SGN X*ABS Хжжү 
Z=truncation to Y plcs.(X) LET Z-SGN X*INT (ABS Xx*x10xxY)?/10xxY 
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Appendix В PHYSICAL AND MATHEMATICAL CONSTANTS 


Values of some of the more or less fundamental constants are 
shown in table B-1. The physical constants are from reference (22). 
Uncertainties are given in parts per million (ppm) of the main 
number. For example, the uncertainty in the value of the speed of 


light is: 


(2. 99792458»10”ኞ)።‹ (.004x10 "ን m/s = +1.2m/s 


‘The mathematical constants are from references (1) and (3). 


::788%06606/529661291977 = ZM 
7776801159/28157<0618%776296%0 = ә ОТЗот 
'''6166:10989650956626058620€'Z = ОТ чт 
::71206041%760657665<081/%766970 = 2 
:7771590909826<106799512/150 
57: 1820966625%066%82819287/:2 
:77%918/06266276%7/60285075/6666697/6778820<6/268Е6%929786766/69С6<926<191:6 


14816000 Тепотівітлело X 02/99 
88вш 16ӘЛ полаоетч 76660176 
29235409 пиеш2а1о9 с7-01 Х 0990851 
30215409 $, моцета uddy*¢ + ус-0Т X 18869601 


еблечо подота 1946 °С + 61-017 * 268120971 


шаповд ит 34877 Jo 9598 801900" + 601 х 85726/6672 


ледшту 8, охреЗолу uddp*c + OT X 47022079 


(sarun IS) твотзАча 


8388938903 1-8 этает, 
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Appendix С UNIT CONVERSIONS 


The following series of conversion charts give conversion 
factors within and between metric and U.S. Customary systems of 
units. All conversion factors are exact except those noted, in 
which case the exact expression is found in the footnote. The 
charts show the units in their relative sizes on a logarithmic 
axis, with the smallest at bottom and the largest at top. For 
example, on the Length chart, there are 100 centimeters (cm) in 
a meter (п), so the label for the cm is spaced 2 units below the 
label for m, since 2 - log, 100. 

To obtain a conversion factor between any two units on a chart: 


1. Locate the two units, UNIT1 and UNIT2,and find a path between them. 


2. Let Е-1 
Move along the path between UNIT1 and UNIT2: 
If moving in direction of arrow, then let F=F/NUMBER ON PATH 
If moving in opposition to arrow, then let F-FxNUMBER ON PATH 
At end of path, check: If UNIT2 is lower than UNIT1, then Е 


should be greater than 1. 

3. Е has units of UNIT2/UNIT1: 
If converting X UNIT1's to Y UNIT2's, then let Y-XxF 
If converting Y UNIT2's to X UNIT1's, then let X-Y/F 


4. Example: Convert 15.23 meters to inches. 


From the Length chart: 


Starting at meters (т) and moving toward inches along the dotted 
path: 


F=1*107/2,54  in/m 


15.23m = 15,23*F іп = 599.61 in 


Conversion Charts notes: 


1. 6336/3937 4. (5659/6640)х1.9 
2. (3937/6336)7x6.4 5. 175/144 
3. .67200625х1.6387064 
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Appendix р SUMMARY ОҒ CENTRAL SUBROUTINES AND PROGRAMS 


Tables D-1 and D-2 1155 the properties of the central sub- 
routines and programs from chapter 1 through chapter 8. The column 
labeled "User-supplied variables" lists the variables which а 
subroutine or program prompts the user for. The column labeled 
"Caller-supplied variables" lists the variables which are supplied 
by a calling main program or subroutine. The column labeled 
"Output variables" lists the variables which a central subroutine 
returns to the calling main program or subroutine. The column 
labeled "Variables changed" lists all variables changed by the 
central program or subroutine, as well as those changed by any 
subroutines which are called by that central subroutine or program. 
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Please use the following page to send in any comments, suggestions 
interesting results, new programs, etc, If you find any glaring 
errors in the book, please correct me. If you have succeeded in 
translating any of these programs to different computers, please 
tell me how you did it, and what problems you had in doing it. 
If there is a second printing, I may include an extra chapter of 
programs sent in by readers. 
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